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ABSTRACT 


This  investigation  explores  the  effects  of  chordwise  forces  and 
deformations  and  steady-state  deformation  due  to  lift  on  the  static 
and  dynamic  aeroelastic  stability  of  a uniform  cantilever  wing.  Results 


of  this  analysis  are  believed  to  have  practical  applications  for  high- 
performance  sailplanes  and  certain  RPV's. 

The  airfoil  cross  section  is  assumed  to  be  symmetric  and  camber 


bending  is  neglected.  Motions  in  vertical  bending,  fore-and-aft 
bending,  and  torsion  are  considered.  A differential  equation  model  is 
developed,  which  includes  the  nonlinear  elastic  bending-torsion  coupling 
that  accompanies  even  moderate  deflections.  A linearized  expansion  in 
small  time-dependent  perturbation  deflections  is  made  about  a steady 
flight  condition.  The  stability  determinant  of  the  linearized  system 
then  contains  coefficients  that  depend  on  steady  displacements.  Loads 
derived  from  two-dimensional  incompressible  aerodynamic  theory  are  used 
to  obtain  the  majority  of  the  results,  but  cases  using  three-dimensional 
subsonic  compressible  theory  are  also  studied. 

The  stability  analysis  is  carried  out  in  terms  of  the  dynamically 
uncoupled  natural  modes  of  vibration  of  the  uniform  cantilever.  Dynamic 
stability  in  the  case  of  incompressible  strip-theory  airloads  is  deter- 
mined in  two  ways.  One  is  the  "V-g  method"  familiar  to  aeroelasticians . 


When  steady  deformations  are  present  this  method  requires  an  iterative 

1 

matching  of  flutter  speeds  with  estimated  speed.  The  second  approach  ' 


• Scctlba 


involves  determination  of  the  complex  eigenvalues  of  the  aeroelastic 


modes  at  any  desired  flight  condition.  The  aerodynamic  loads  are 
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expressed  in  terms  of  the  generalized  Theodorsen  function;  eigenvalues 
of  the  aeroelastic  system  are  located  with  a gradient  search  technique. 

The  effect  of  steady  drag  on  flutter  of  nonlifting  wings  using 
incompressible  strip-theory  is  studied  and  shown  to  correlate  with 
previously  known  results.  Next,  the  influence  of  steady  lifting 
deformations  on  flutter  is  investigated,  and  flutter  modes  are  found 
that  involve  fore-and-aft  bending  motions.  The  significance  of  unsteady 
leading  edge  suction  forces,  which  are  predicted  by  the  two-dimensional 
incompressible  aerodynamic  theory,  is  then  examined.  Two  idealized 
examples  based  upon  existing  sailplanes  are  analyzed. 

Steady  drag  loads  lower  the  flutter  speed  for  larger  aspect  ratios 
but  increase  it  for  aspect  ratios  less  than  a certain  value.  Divergence 
speed  is  more  sensitive  to  steady  drag,  and  for  very  high  aspect  ratio 
wings  it  can  fall  below  the  bending-torsion  flutter  speed.  Steady 
deformations  due  to  lift  always  decrease  the  flutter  speed  by  an  amount 
dependent  upon  the  aspect  ratio  and  the  fore-and-aft  bending  stiffness. 
Leading-edge  suction  forces  increase  flutter  speed. 

Three-dimensional  steady  and  unsteady  airloads  are  introduced  into 
the  V-g  flutter  analysis  scheme,  and  for  a Mach  number  of  zero  the  role 
of  steady  lifting  deformations  and  unsteady  leading-edge  suction  forces 
is  more  accurately  determined.  The  behavior  predicted  using  strip 
theory  loads  is  again  observed,  and  the  suction  forces  are  confirmed  to 
contribute  a significant  stabilizing  effect.  Further  calculations  using 
high  subsonic  Mach  numbers  reveal  only  mild  effects  due  to  compressi- 
bility (disregarding  unsteady  chordwise  loads  of  viscous  origin). 
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Chapter  I 
INTRODUCTION 

Aerodynamic  loads  and  deformations  parallel  to  the  chord  are 
usually  neglected  during  studies  of  aeroelastic  stability  of  lifting 
surfaces.  Furthermore,  dynamic  stability  is  usually  analyzed  without 
regard  for  the  influence  of  steady  deformations  due  to  steady-state 
lift. 

Very  little  literature  exists  which  treats  chordwise  forces  and 
bending  in  an  aeroelastic  analysis.  A substantial  study  of  the  effects 
of  drag  loads  on  divergence  of  a cantilever  wing  is  made  by  Petre  (Ref. 
2,  pp.  449-487).  Here  it  is  clearly  demonstrated  that  the  interaction 
of  drag  with  bending  deformations  due  to  lifting  loads  can  signifi- 
cantly reduce  divergence  speeds.  Goetz  (Ref.  17)  considered  this  same 
drag-bending  deformation  divergence  mechanism,  specialized  to  the  case 
of  a rigid  lifting  surface  at  the  end  of  a beam-rod.  This  work  involved 
supersonic  flow  past  a surface  having  a blunt  leading  edge,  and  the 
resulting  sizeable  drag  forces  caused  a significant  reduction  of  the 
classical  divergence  speed. 

One  example  of  an  aeroelastic  study  in  which  chordwise  deformations 
of  a straight  cantilever  wing  are  accounted  for  is  the  work  in  the  area 
of  tilting  proprotor  aircraft  by  Wayne  Johnson.  The  cruising  flight 
condition  (Ref.  27)  is  modeled  using  a proprotor  with  axial  flow  mounted 
at  the  tip  of  a cantilever  wing.  The  additional  degrees  of  freedom 
associated  with  the  individual  elastic  rotor  blades  and  the  aerodynamic 
and  inertial  effects  of  the  proprotor  result  in  a much  more  complex  and 
specialized  aeroelastic  analysis  than  is  considered  here. 
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The  effect  which  steady  chordwise  forces  can  have  upon  dynamic 
stability  was  explored  by  Petre  (Ref.  3)  and  Petre  and  Ashley  (Ref.  1) 
using  two-dimensional  incompressible  unsteady  lifting  airloads.  The 
latter  work  presents  extensive  calculations  regarding  the  effect  of 
steady  drag  on  bending-torsion  flutter  of  a uniform  cantilever  wing. 

It  serves  as  the  starting  point  for  the  work  pursued  in  this  thesis. 

The  objectives  of  this  thesis  are  as  follows: 

1.  To  check  and  interpret  the  predicted  effect  of  steady  drag 
on  the  flutter  behavior  of  a nonlifting  wing  discovered  in 
Ref.  1,  using  a modal  approach  instead  of  a collocation 
approach . 

2.  To  generalize  the  equations  of  motion  to  include  consistently 
fore-and-aft  bending  motions,  adequately  accounting  for  the 
elastic  coupling  among  the  three  degrees  of  freedom. 

3.  To  include  steady-state  lifting  deformations  in  the  dynamic 
stability  analysis  by  considering  small  time-dependent  per- 
turbation deflections  about  a steady  displacement  solution. 

A.  To  allow  for  unsteady  leading-edge-suction  forces  in  the 

chordwise  direction  predicted  by  two-dimensional  incompress ibl e 
unsteady  potential  flow  theory. 

5.  To  improve  the  representation  of  both  steady  and  unsteady  air- 
loads by  use  of  a three-dimensional  subsonic  kernel  function 
program,  from  which  leading-edge  suction  and  induced  drag  can 
also  be  obtained. 

Items  2,  3,  and  A are  interrelated  and  together  represent  a consis- 
tent extension  over  previous  research  in  the  modeling  of  the  physical 
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system.  For  example,  when  steady  lateral  bending  deformations  are 
present,  unsteady  chordwise  loads  can  induce  twisting  motions  which 
significantly  affect  dynamic  stability.  The  main  purpose  here  will 
be  to  establish  trends  and  gain  fundamental  insight  into  the  influence 
of  chordwise  forces  and  steady  deformations,  hopefully  shedding  light 
on  their  importance  in  practical  aerospace  problems. 

Certain  assumptions  are  adhered  to  throughout  this  thesis.  Since 
the  emphasis  is  upon  working  from  the  equations  of  motion  in  differen- 
tial form  in  order  to  include  certain  nonlinear  elastic  coupling  terms, 
it  is  convenient  to  restrict  this  study  to  straight  cantilever  wings 
having  mass  and  stiffness  properties  uniform  with  span.  The  wing  is 
taken  to  be  a one-dimensional  structure  in  the  sense  that  all  deforma- 
tions are  described  as  functions  of  the  spanwise  variable  y . Camber 
bending  is  neglected  and  the  simple  Euler-Bernoulli  beam  stress-strain 
assumptions  are  used.  The  platform  is  rectangular,  and  the  steady  and 
unsteady  flow  fields  are  always  assumed  to  be  superposable;  unsteady 
loads  are  computed  for  the  undeformed  geometry  and  applied  to  the  deformed 
wing.  Although  these  assumptions  would  be  restrictive  for  the  purpose 
of  modeling  actual  structures,  they  are  acceptable  here  since  only  the 
relative  influences  of  chordwise  loads  and  steady  deformations  are  of 
interest . 

In  Chapter  II  a vertical-bending/torsion  model  basic  to  the  system 
of  Ref.  1 is  developed  and  modal  equations  are  derived  to  permit  flutter 
calculations  for  zero  steady  lift  with  steady  drag  included.  Assumed 
mode  solutions  are  then  compared  with  results  of  Ref.  1.  A linear 
steady-state  version  of  the  modal  equations  is  then  examined  to  allow 
computation  of  divergence  speeds  as  affected  bv  steady  drag. 
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In  Chapter  III  a model  central  to  this  thesis  is  developed  that 


includes  vertical  bending,  chordwise  bending,  and  torsion.  Nonlinear 
steady  and  linearized  unsteady  differential  equations  are  then  deduced, 
and  they  are  analytically  compared  with  the  model  in  Chapter  II. 

In  Chapter  IV  the  modal  forms  of  these  steady  and  linearized 
unsteady  equations  are  set  up  to  include  lifting  airloads  deduced  from 
incompressible  steady  and  unsteady  strip-theory.  A scheme  based  on  the 
so-called  V-g  method  of  flutter  analysis  is  used  to  determine  neutral 
dynamic  stability  conditions,  and  results  are  checked  against  those  of 
Chapter  II. 

In  Chapter  V'  a generalization  to  the  case  of  arbitrary  motion  in 
time  is  presented  through  Laplace  transformation  of  the  modal  equations, 
which  requires  incompressible  unsteady  two-dimensional  airloads  valid 
for  non-periodic  motions  of  the  wing. 

A determinant  iteration  procedure  is  used  to  determine  the  aero- 
elastic  eigenvalues  for  flight  speeds  above  and  below  the  flutter  speed. 
Finally,  the  effect  of  unsteady  leading-edge  suction  forces  predicted 
by  incompressible  strip  theory  is  included  in  the  linearized  unsteady 
stability  system. 

In  Chapter  VI  all  results  for  incompressible  strip-theory  airloads 

✓ 

are  assembled  and  systematically  presented,  concluding  with  two  examples 
based  upon  actual  sailplanes. 

In  Chapter  VII  the  flutter  speed  prediction  scheme  of  Chapter  IV 
is  modified  to  use  three-dimensional  subsonic  steady  and  unsteady 
lifting  airloads.  Results  are  presented  to  indicate  the  effects  of 
three-dimensional  aerodynamics,  unsteady  drag,  and  compressibility. 
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Chapter  II 

UK  AC  EFFECTS  ON  FLUTTER  AND  DIVERGENCE; 

DEFORMATION  !N  VERTICA1.  BEND INC  AND  TORSION  ONLY 

A.  The  Vertical  Bend  i ng/Tors  i on  Equations 

In  the  Petre  and  Ashley  work  (Ref.  1)  the  effect  of  drag  on  flutter 
of  cantilever  wings  Is  studied  assuming  only  vertical  bending  and 
torsion  about  the  elastic  axis.  The  equations  are  cast  into  integral 
form,  and  solutions  are  obtained  by  collocation  of  the  integral  equations 
at  ton  stations  across  the  span.  Here  we  check  the  results  of  Ref.  1 
using  an  entirely  different,  modal  approach.  The  results  also  serve  as 
examples  with  which  to  compare  solutions  found  with  the  more  general 
system,  developed  in  the  following  chapter,  which  includes  ehordwiso 
bend  I ng. 

The  differential  equations  of  bending  and  torsion,  as  given  in 
Ref.  1 and  adapted  to  the  present  notation  and  coordinate  system  (as 
In  Fig.  2-1),  are 


(2- la)  | F.I  -ir—rl  + m 7 - s ly  4-  - 1.  (w,*f>) 

av  x Ay  at  e at  a 


+ -K-7  (M  <f>)  - 0 
Av  7. 


(2-lb) 


§7  "Ti  $1 


at-  ' ”0  at'  ’ ma(w,<^) 


3' w 

- M -^-9  - 0 

7.  av 


7 


FAQ* 

\ ' 


I 


The  quantity  is  a function  of  y given  by 

(2-2)  M - - / D(n) (n-y)dn 
z y 

= - d st^-y 1 2 

where  it  is  assumed  that  the  running  drag  force  D has  constant  magni- 
tude across  the  span.  M can  be  recognized  as  the  total  moment  about 

z 

the  vertical  axis  applied  at  station  y by  all  drag  acting  outboard 
of  this  station.  Mean  values  of  the  twist  and  bending  displacement 
are  assumed  to  be  zero. 

Equations  (2-1)  are  linear  and  are  coupled  inertially  by  the  third 
term  in  each;  these  terms  arise  from  the  offset  of  the  center  of  mass 
of  the  airfoil  section  from  the  elastic  axis.  A second  coupling  effect 
is  due  entirely  to  drag,  introduced  through  the  terms  containing  M 

z 

All  remaining  terms  in  the  equations  represent  the  conventional  elastic, 
inertial,  and  unsteady  lifting  aerodynamic  load  contributions. 

The  manner  in  which  the  drag  loads  have  been  introduced  into  the 
system  is  discussed  by  i’etre  (pp.  449-487,  Ref.  2)  and  can  be  explained 
by  the  following  physically  oriented  argument.  The  drag  coupling  term 
in  the  bending  equation  arises  from  the  resultant  bending  moment  m^(y) 
applied  at  station  y duo  to  drag  forces  outboard  of  v . As  shown  in 
Eig.  2-1,  a drag  force  acting  outboard  of  station  v has  a component 
l*Mv)  perpendicular  to  the  airfoil’s  principal  axis  of  vertical  bending, 
giving  a resultant  moment  at  y about  the  principal  direction  of 

m (y)  - <Ky)n  / (n-y)dn 
x y 

- <fr(y)D  S(e-y)2 

- -M  4> 
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The  moment-curvature  relation  for  the  beam. 
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together  with  moment-shear  equilibrium 
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9 7s"  = q(y) 

(here  q(y)  is  the  positive-upward  running  load  in  the  principal 
direction  of  the  section)  allows  the  drag  effect  to  be  expressed  in 
equilibrium  with  the  elastic  term  as 


9 2 
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This  is  the  same  drag  coupling  term  appearing  in  (2-1). 

A similar  derivation  reveals  the  origin  of  the  drag  term  in  the 
torsion  equation.  Looking  at  the  front  view  of  the  wing  (Fig.  2.2)  one 
can  see  that  a drag  load  at  n gives  rise  to  a twisting  moment  at  y 
by  acting  through  the  moment  arm  given  by  distance  e , 


e = w(n)  - £w(y)  + ^ (y)(n-v)j 


The  resultant  torque  applied  at  station  v by  all  drag  forces  acting 
outboard  is  given  by 


t (y)  = d [w(n)  - w ( v ) - ^V')  (n-v)  Idn 

y y <Jy 


Differentiating  with  respect  to  y . 
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Elastic  equilibrium  for  a rod  loaded  with  applied  torsion  Ty(v)  is 
given  by 
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The  drag  coupling  term  of  equation  (2-1)  can  then  be  identified  in 
the  result 


9_ 

9y 


= 0 


B.  Solution  By  Assumed  Modes 

When  structural  dynamics  problems  yield  solutions  whose  frequencies 
are  within  the  range  of  the  structure's  lowest  normal  mode  natural 
frequencies  and  the  latter  have  a sufficiently  sparse  distribution, 
modal  analysis  methods  prove  to  be  effective.  Primary  bending-torsion 
flutter  of  cantilever  wings,  under  study  here,  is  a classic  example  in 
aeroelasticity  of  such  a system.  Inclusion  of  steady  drag  effects 
should  have  only  a small  effect  on  the  range  of  frequencies  over  which 
flutter  solutions  occur  and,  by  this  reasoning,  should  not  adversely 
affect  the  convergence  of  modal  solutions.  Petre  in  Ref.  1 expresses 
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the  opinion  that  methods  relying  on  the  assumption  of  a few  prescribed 
deformation  modes  constrain  the  flutter  mode  shape  and  should  be  avoided. 
This  opinion  is  tested  herein  by  actually  studying  modal  convergence. 

Solutions  of  (2-1)  are  sought  using  superposition  of  the  dynami- 
cally uncoupled  natural  modes  of  free  vibration  in  bending  and  torsion 
of  a uniform  cantilever  beam.  These  assumed  modes,  with  their  relevant 
properties,  are  described  in  Appendix  A.  Although  not  true  normal  modes 
of  the  inertially-coupled  structure,  they  can  be  considered  as  "pseudo- 
orthogonal"  since  integrals  of  the  type 


Ji  m f f dv  = 0 
0 wi  w. 

;oJ  \ ^ = ° 


1 * .1 

i * .1 


lead  to  uncoupled  elastic  behavior  and  hence  a diagonal  stiffness  matrix 
in  the  matrix  eigenvalue  problem.  Use  of  the  actual  normal  modes  would 
require  that  they  be  calculated  for  each  wing  configuration  studied. 

Since  the  assumed  modes  satisfy  the  natural  boundary  conditions  at  the 
free  end  of  the  cantilever  as  well  as  the  geometric  boundarv  conditions 
at  its  clamped  root,  C.alerkin's  method  can  be  applied  to  the  differential 
equations  (Ref.  4,  p.  218)  to  obtain  the  system  in  terms  of  modal 
generalized  coordinates. 

To  find  the  velocity  for  neutral  stability  (the  flutter  velocity), 
the  V-g  method  (p.  381,  Ref.  5),  common  in  aeroelastic  stability 
analyses,  in  employed.  With  simple  harmonic  motion  of  frequency  u) 
the  unsteady,  incompressible,  strip-theory  lifting  airloads  are 
expressible  as  [(4-123)  and  (4-124)  of  Ref.  51 
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(2-3a)  L (wj;t)  = — IT  p b 3 io2{-  L £ + [L.  -('s+a)L  ]4>}e1(i)t 
d w b <p  w 

(2-3b)  m (w,4>;  t)  = II  p b4  u)2{-[M  -(!j+a)L  ] £ 

cl  W w D 

+ [M^  - (Jjfa)(L^  + hw)  + (h+a)2  Lw]^}ei(l>t 

where 

/L' 

(2-4) 

L 

( 
k 

Here  C(k)  is  the  familiar  Theodorsen  function  of  reduced  frequency 
. cob 

k " T • 

The  drag  coupling  terms  in  (2-1)  can  be  treated  as  applied  loads 
in  developing  the  modal  equations,  by  defining  total  applied  force  and 
moment  in  the  bending  and  torsion  equations  as 


(2-5a) 

F (y » t)  = L (w,<J),  t) 

^ cl 

■ fy*  (Mz(J>) 

(2- 5b) 

my(y,t)  = ma(w, if),  t) 

32w 

-m 

z dy 

■ 1 - T c(k)  \ - H 

■ H-i  [1  +2C(k)J  -^C(k) 


Incorporation  of  artificial  structural  damping  g by  allowing  a complex 
elastic  modulus  produces  the  system 


(2-6n)  Eyi+1,)  + » Ip  - 0 ■ yy.t) 

(2-1  ) OId(l+ig)  + my(y,t) 
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Formal  development  of  the  modal  equations  begins  with  substitution 
of  the  series  expansions  for  w and  <f>  in  terms  of  the  assumed  mode 
shape  functions  and  generalized  coordinates.  With  the  same  number  n 
of  bending  and  torsion  modes  always  used,  the  system  order  will  be  2n  . 
Generalized  displacements  are  assumed  to  have  simple  harmonic  time 
dependence,  giving 


(2-7a)  w(y , t)  = £ f (y)  q e 

i=l  i Wi 


iwt 


(2- 7b)  <t>(v,t)  = E f (v)  q e 

1=1  *i  ' *i 
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Generalized  displacements  for  bending  modes  have  units  of  length,  where 
as  for  torsion  they  are  dimensionless. 

Galerkin's  method  involves  substitution  of  (2-7)  into  the  system 
represented  by  (2-1),  (2-2),  (2-3),  and  (2-4)  and  then  multiplication 


of  each  term  in  the  bending  equation  by  fw^  and  each  in  the  torsion 


equation  by  f^  , followed  by  integration  across  the  half  spa 


alf  span.  With 


iuit 


e cancelled  the  resulting  system  is 
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(1  £ j £ n) 
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The  generalized  forces  can  be  arranged  in  terms  of  dimensionless  unit 

generalized  forces  Q through  the  following  definitions: 

r cj 

(2-9a)  A (<J>,w;t)f  dy  = ITpu)V)J  Z Q -r1  + E Q,  ..  q.  ela 
02  wj  [_1=1  J’1  b i=i  J*1+n  ^ 


w,  n 


(2-9b)  jTm  (*,w;t)f . dy  = IIpwVl  T,  Q , + E Q . . q. 

0 y 1=1  j+n , i b 1=1  j+n,i+n  4^ 


The  modal  integral  in  the  first  term  of  (2-8a),  integrated  by  parts 
twice  and  with  application  of  cantilever  boundary  conditions,  yields 


/ f""  f dy  = / f”  f"  dv 
o wt  w|  0 w^  w^ 


0 for  i ^ j 


From  Appendix  A,  the  i=j  term  can  be  expressed  in  terms  of  the  natural 
frequency  of  the  jth  assumed  bending  mode  by  the  substitution 

(2-10)  El  jT*(f")2dy  = mu2  f2  dy  = mu>2  l 
x 0 w n ... 


W . 0 w^ 


Similarly,  integrating  the  first  term  in  (2-8b)  by  parts  and  introducing 
the  natural  frequency  (0^  of  the  jth  assumed  torsion  mode  leads  to 

(2-11)  Cl.  / f"  f.  dy  = f 

a o <J>i  <J>^  / 0 (i  r j) 

- ?s  Jto2  l (i  = j) 


Insertion  of  (2-9),  (2-10),  and  (2-11)  into  (2-8)  and  further  use 
of  modal  orthogonality  properties  gives 
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(1  _<  j _<  n) 
(2-12b) 


t 1a  < 

1 m J d <p 

2 JTpb2  mb"  J to2 


s n w. 

, _e  m , t i 

bm  npb2  ^ J j ij  b 


- U. 


npb2  mb2  q<j 


.1 


+ 1 Q 

i = I 


.i+n,i  b 


+ T.  Q. 


1 = 1 i+n«1+n  -i-j 


<t>,  = 0 


where  the  bending  equation  has  been  divided  through  by  ITpu)2b3t  and 
the  torsion  equation  by  IIpu^bH  . Inertial  coupling,  a consequence 
of  nonorthogonality  between  bending  and  torsion  assumed  modes,  produces 
terms  in  the  modal  integrals 


(2-13)  I = fl  f (v)f.  (v)dv 

i i o w.  • (p . 

i j 


For  now  it  is  convenient  to  reference  to  to  the  first  assumed 


torsion  mode  natural  frequency  to 


Assumed  mode  natural  frequencies 


(Appendix  A)  are 


to 


(2-14) 


2 i n2  “Ld 

4-j'  2 V j 


in'* 

,2  »T  4 IT  X 

to  N , - 

w . i It  m 

i 


Ratios  that  are  useful  are  then 
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(2-15) 


Cl)  1.  0 

W,  N^n-  EI  b2  _ 

_Li  = _j_  x j 

< 4 gTTF  mb' 

0,  d 


IT1  = C2  j-l) 


The  transcendental  are  given  in  Appendix  A for  1 < j < 5 . 

The  equations  can  be  nondimensionalized  through  the  following 
dimensionless  parameters,  chosen  to  be  as  consistent  as  possible  with 


Ref.  1 notation. 


(2-16) 


s 


a ~ mb' 


Frequency  is  nondimensionalized  as  in  Ref.  1 by  defining  the  frequency 


parameter 


(2-17)  Z = (l+ig)ft2  = (1+ig)- 

u)2 

4 *1 


7F  Td7-  (1+ig) 


Substitution  of  (2-15),  (2-16),  and  (2-17)  into  (2-12)  yields  the  modal 
equations  in  a suitable  form  for  computation. 

q q 

w.  n Mw. 

(2-18a)  (M  - Nfn-MPi  Zh-r1  + E Q — - 

1 a b j.j  J,i  b 


(1  < j < n) 


+ £ {MS’jl  + ‘WS,  ’ ° 
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(2-18b) 


w 


MS,ii  + V.i1  IT  + ^ V1  - n’O-W'zi)^ 


+ T.  Q. 


1=1  '‘j-Hi.i+n  q4>1  = ° 


The  flutter  determinant,  of  order  2n  , results  as  a necessary  condition 
for  a nontrivial  solution  and  leads  to  a complex  eigenvalue  problem  to 
determine  Z . 


Next,  the  unit  generalized  forces  . must  be  expressed  to  allow 


numerical  computation.  Combining  (2-3),  (2-4),  and  (2-7)  into  (2-9) 
leads  to 

% 

i 

: — dv  + . 

0 


(2-19a)  - IlpbV  {-/o£  Lw  E fw  (y)fw  (y)  ~ dv  + f}  [LA  - (V a)Lj 

i=l  i .i 


i:i%i(y)fw.(y)Vy}  + fo  ^ [l  V2bCD(t-y)2iSif(|(i(y)q^: 


(1  < j < n) 


• f (y)dv  = npu»2b*4[  I Q . + Z Q.  . . q.  ] 

| i=1  J.1  b 1=1  .1,i+n  ^<t> 


w^  n 


n n W. 

(2-19b)  lTpbVf  r [-M  + Ci+a)L  ] Z f (v)f  (v)  -r1  dv 

0 w w wi  ^ 


+ fo[n6  ~ (VHa)(L  +Mw)  + (^a)2Lw]  Z f (y)f  (y)q  dy> 


+ f}  h PV21 


i 


f0  dy 


= npto2b4e(  e o 


w n 

+ E 0. 


1 = 1 I411*4  b i = l ,1+n.l+n 


After  further  reduction,  with 


,2  u>‘b‘ 
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I 


and  carrying  out  remaining  differentiations  in  y , (2-19)  becomes: 


(2-2°a)  (y)dy)  -p  + + 'H+.KJ  • 


• f„/>d*  + !ipi2  '\$K  <«d?  - 4 • 


(1  < j < n) 


• f!  (y)fw  (y)dy  + £ d-y)2f"  Cy)fw  (y)dy]}q 

i j i j Ti 


n w . n 

E Q,  . “7“  + E Q.  . q, 
i_l  J*1  b i=1  j , l+n  >i 


(2-20b) 


E {_MW  + (^+a)V^lfw^(y)f<^(y)dv 


i=l 


+ w /oL(1-y)2fw/?)f<}».(y)d?}~bi  + <E1{[v(,'5+a)(VMw) 

i J i=l 


+ (!’+a)  2L  ] f}f.  (y)f,  (y)dv }q 


w 1 0 <J>.  J q>. 

1 j 


♦i 


n 


n 

if1  Qj+n,i  b + iE1  ^j+n,i+n  q<t>1 


Here,  integrals  in  the  spanwise  variable  y have  been  non- 
dimensionalized . Four  different  types  of  modal  integrals  are  encountered, 
including  the  1^  previously  identified.  The  three  new  forms  arise 
from  the  drag  coupling  terms  and  are 

/ 


I 


if  ’ /01(1‘5)2fW1<>i)f».<5)d? 


(2-21)  < 


’if  - 


I' 


V -if  ■ (,1<l-v)2f;i<y)fUj«)0v 
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The  modal  integrals  in  (2-13)  and  (2-21)  were  numerically  computed  to 
twelve  significant  digits  for  1 _<  i,j  £ 5 . 

Two  additional  dimensionless  parameters  can  be  introduced,  again 
drawn  from  Ref.  1. 

C 

(2-22)  A = h + 2 C = W 

The  drag  parameter  C is  defined  as  the  applied  steady  sectional  drag 
coefficient  divided  by  the  sectional  lift  curve  slope  (211  for 
incompresible  strip  theory). 

After  (2-20),  (2-21),  and  (2-22)  are  combined  with  orthogonal 
modal  integrals  recognized,  the  unit  generalized  forces  are  found  to  be 


(2-23) 


> 

= 1 L 

(i=i) 

< W 

1° 

(i*j) 

i 

j , i+n 

= ~ [L 

- AL  ]I 
w 

^j+n,i 


F [2lii 


[M  - AL  ]I. . + ~r  lf1} 
w wJ  i,i  k l j 


4 I(2)  + 1(3) 
i.l  tJ 


^j+n, i+n 


^ + V + A V (i=J> 

0 («j) 


The  Qj.  depend  only  upon  k.  A,  and  C . Actual  computation  of  the 
Theodors en  function  is  accomplished  by  direct  use  of  the  ascending 
power  series  of  the  modified  Bessel  functions  and  , as  explained 

in  Appendix  B. 

Equations  (2-4),  (2-18),  and  (2-23)  together  supply  the  flutter 
determinant  for  the  VT  system  using  assumed  modes,  given  in  matrix 
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I 

form  in  (2-24).  Solutions  can  be  found  using  from  one  to  five  assumed 
modes  in  each  of  bending  and  torsion,  and  the  maximum  order  of  the 
determinant  will  be  10. 


(2-24) 


M[1  - + qu  T0~l 

L°-^M[l  -(N  n)4Pi  Z]  + Q 

n a 'n , n 


i 

1 -MSI, , +Q 
i 11  'l,n+l 

I 

I 

i ! 

1 -MSI  +Q 
j ni  n,n+l 

I 


MSI 


+ 0, 


In 


1, 2n 

i 

i 

-MSI 

nn 

+ (^n,2n 


-MSI  + Q , . . 
11  n+1,1 


MSI  . + Q . , 1 hMi  [1-  "Z]  + 0 ^n"1 

nl  n+1 ,n  | a 4 n+1, n+1 


■MSI  + 0,  . 

In  2n,l 


-MSI  + Q9 

nn  2n,n 


2tt2  , 


In  seeking  solutions,  values  for  M,  P,  i^.  A,  and  S must  first 
be  chosen  to  specify  the  wing  configuration.  Then,  for  any  desired 
value  of  C and  an  estimate  of  reduced  frequency  k , complex  eigen- 
values Z can  be  found  from  (2-24)  via  linear  matrix  eigenvalue  techniques. 

A computer  program  is  used  to  solve  for  flutter  conditions  as 
follows.  With  a first  estimate  of  k chosen  large  enough  so  that  the 
structural  damping, 


(2-25) 


Im  (Z) 
8 = Re(  Z) 


is  negative  for  all  2n  eigenvalues,  successively  smaller  values  of  k 
are  assigned  and  eigenvalues  computed  until  a positive  g is  obtained 
for  the  eigenvalue  corresponding  to  the  aeroelastic  mode  which  encounters 
flutter.  Then  a zero-finding  subroutine  locates  k for  which  g=0 
for  the  flutter  mode.  Dimensionless  speed  and  frequency,  defined  by 
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('.  Compar  isoji  of  Mocl.il  with  Co  I I oo.it  1 on  Results 

A rather  thorough  study  was  carried  out  In  Ref.  1,  covering  a wide 
range  of  practical  combinations  of  (lie  dimensionless  parameters.  The 
present  objective  Is  not  to  recalculate  all  of  the  same  data  hut  rather 
to  evaluate  the  effectiveness  of  the'  assumed-modi'  approach.  Consequently, 
(2-24)  has  been  solved  at  conditions  parallel  to  ones  for  which  results 
are  published  In  Ref.  1 to  offer  direct  compari son . Also  (2-24)  Is  used 
to  verify  the  performance  of  the  VCT  system  developed  in  the 
next  chapter. 

A comparison  of  assumed  mode  calculations  with  Ref.  1 results  is 

presented  in  Table  2.1  for  three  different  configurations.  Of  these, 

cases  (a)  and  (c)  represent  a stubby  low-aspect-ratio  wing  (with  T ” .4) 

whereas  case  (b)  Is  the  opposite  extreme  of  large  aspect  ratio.  A 
El 

x 

typical  ratio  -yr—  = 1.6  , for  example,  would  fix  the  aspect  ratio  of 
il 
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cases  (a)  and  (c)  at  2 and  of  case  (b)  at  20  . The  mass  ratio  para- 
meter of  cases  (a) , (b)  is  in  the  regime  encountered  by  sailplanes 
whereas  in  (c)  it  is  representative  of  conventional  aircraft.  Finally, 
the  steady  drag  parameter  is  zero  in  the  first  two  cases  but  has  the 
extremely  large  value  C = 0.04  in  (c)  . 

In  all  cases  n=3  already  yields  adequate  modal  convergence. 

For  the  large  aspect  ratio  case  the  use  of  one  assumed  mode  (n=l) 
gives  a significant  error,  yet  n=2  produces  good  accuracy.  This 
suggests  that  the  second  bending  mode  is  an  important  factor  in  flutter 
of  high-aspect-ratio  wings,  a phenomenon  discussed  in  Chapter  b. 

Additional  comparisons  between  flutter  speeds  and  frequencies 
found  by  the  two  methods  are  offered  in  Table  2.2,  emphasizing  their 
relative  accuracy  as  steady  drag  is  increased  to  the  very  high  value 
C = 0.04  . Evidently  good  agreement  is  maintained  in  the  presenci  >f 
drag. 


The  mode  shape  at  flutter  for  two  of  the  preceding  cases  is 

presented  in  Table  2.3  for  wings  of  small  and  large  aspect  ratio,  each 

with  C = 0.04  . Amplitudes  of  the  ten  generalized  displacements  are 

normalized  with  respect  to  | q | , and  the  phase  angles  are  referenced 

4>1 

to  the  phase  of  this  torsion  mode.  For  F = 0.4  only  q.  and  q 

*1  W1 


have  appreciable  magnitude,  whereas  the  large-aspect-ratio  example 

displays  a significant  contribution  bv  q as  well. 

W2 

These  flutter  mode  shapes  can  be  compared  to  the  Ref.  1 results, 
with  some  effort,  as  follows.  In  Ref.  1,  the  tangent  of  the  angle  by 
which  the  torsional  displacement  at  the  wingtip  leads  bending  displace- 
ment in  the  flutter  mode  is  tabulated;  it  is  converted  here  to  an  angle 


j 
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in  degrees.  The  spanwise  shape  of  the  bending  portion  of  the  mode  at 
flutter  for  the  two  Table  2.3  cases  is  illustrated  in  Fig.  5 of  Ref.  1. 
Since  the  flutter  mode's  bending  displacement  obtained  in  that  treatment 
was  complex,  its  amplitude  at  each  spanwise  station  was  taken,  phase 
differences  being  neglected  in  the  figure. 

For  the  present  solutions  a similar  assumption  is  made  to  display 
the  bending  mode  shape  of  the  P = 0.004  case  in  Table  2.3.  Contribu- 
tions of  q and  q have  been  added  vectorically  at  points  along 
W1  W2 

the  span  to  allow  comparison  of  bending  amplitudes  with  the  Ref.  1 
figure.  For  p = 0.4  , of  course,  only  the  first  bending  mode  contributes 
significantly  and  phase  differences  are  negligible. 

The  comparison  of  bending  flutter  mode  shapes  appears  in  Fig.  2-3, 
where  amplitudes  are  normalized  to  unit  torsional  displacement  at  the 
tip.  Phase  angles  6 between  bending  and  torsion  at  the  tip  are  also 
compared.  At  low  aspect  ratio  excellent  agreement  for  phase  angles  and 
mode  shapes  is  observed,  with  mild  disagreement  in  the  bending  mode 
amplitudes.  In  the  P = 0.004  case,  for  which  the  Ref.  1 solution  was 
made  using  only  five  spanwise  collocation  points,  a significant  disagree- 
ment in  mode  shapes  and  tip  phase  angles  is  observed.  In  spite  of  this 
discrepancy,  the  respective  flutter  speeds  differ  by  only  0.82%.  For 
all  large-aspect-ratio  cases  compared  this  sort  of  discrepancy  in  flutter 
mode  shapes  is  observed,  where  the  second  bending  mode  plays  a significant 
role  and  causes  appreciable  phase  differences  in  bending  deflections 
along  the  span.  Mode  shapes  at  flutter  are  analyzed  in  more  detail 
in  Chapter  6. 
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D.  Divergence  Including  the  Effects  of  Drag 

One  further  application  of  the  VT  linear  model  is  predicting 


the  effect  of  drag  forces,  steady  in  direction  and  magnitude  and  uniform 
across  the  span,  on  static  aeroelastic  stability.  Removal  of  all  time- 
dependent  terms  in  (2-1)  and  insertion  of  simple  strip-theory  steady 
incompressible  aerodynamic  loads  gives 


o*  w 


(2-28a) 


(2-28b) 


El 


x 9yH 

32<J 


2I7pV2b4>  + i— 7 (M  <t>  ) = 0 

o dy  z o 


32 


GI 


. ~ + 2IIpV‘VA<i)  - M . 

d dy  o z dy 


w 

7^=0 


Subscripts  emphasize  that  deflections  are  static  quantities. 

As  in  the  foregoing  dynamic  analysis,  deflections  are  represented 
by  assumed  modes  and  a system  of  2n  modal  equations  is  derived,  which 
has  a nontrivial  solution  only  if  the  determinant  of  the  matrix  of 
coefficients  is  zero  which  yields  the  divergence  speed  with  drag  effects 
included.  For  brevity,  since  the  manipulations  involved  are  quite 
straightforward,  the  final  form  of  the  stability  determinant  is 
presented  here. 
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Solution  of  this  determinant  for  its  largest  eigenvalue  gives  the 
divergence  speed  . This  is  of  interest  in  the  ensuing  work, 

particularly  for  studying  dynamic  stability  of  wings  having  steady- 
state  lifting  deflections  in  the  presence  of  drag,  which  can  reduce 
divergence  speed  considerably. 

For  zero  steady  drag  (C=0)  the  classical  divergence  speed  comes 
trivially  from  (2-29),  which  degenerates  for  C = C . The  classical 
divergence  mode  shape  is  just  the  first  assumed  torsion  half-sine  mode, 
and  the  first  torsion  modal  equation  uncouples  to  give  the  divergence 
speed, 

(2-30)  UD  = n 

Wien  A is  decreased  to  zero  this  classical  divergence  speed  becomes 
infinite;  yet  when  C / 0 (2-29)  will  still  yield  finite  solutions 

with  A = 0 . 

A direct  comparison  of  solutions  of  (2-29)  with  Ref.  1 results  for 
the  effect  of  drag  on  divergence  is  presented  in  Table  2. A.  Modal 
convergence  is  satisfactory,  but  the  modal  approach  appears  to  differ 
more  significantly  from  the  Ref.  1 analysis  for  divergence  calculations 


than  for  flutter  results. 


Case 


n 


M = 10.0 


a 

S = 0.1 
A = 0.1 
C = 0.0 


TABLE  2.1 


(Ref.  1) 


(Ref.  1) 


(Ref.  1) 


UF 

^F 

% Difference 
of  from  n=5 

2.7175179 

2.7239548 

2.7240004 

1.3105289 

1.3114559 

1.3114641 

0.239% 

0.00239% 

0.00072% 

2.7240178 

2.7240199 

1. 3114673 
1.3114675 

0.000077% 

2.699 

1.309 

0.918% 

4.2621908 

4.0842768 

4.0864182 

4.0866066 

0.842707 

0.8849367 

0.8850560 

0.8850659 

4.296% 

0.0576% 

0.0052% 

0.000597% 

4.0866310 

0.8850660 

— 

4.032 

0.886 

1.337% 

4.260823 

4.260879 

4.260882 

1.294037 

1.2940232 

1.2940250 

0.00155% 

0.000235% 

0.00016% 

4.260889 

1.2940236 

— 

4.199 

1.293 

1.4525% 

Comparison  of  Assumed  Mode  and  Collaeation  Methods  for 
Predicting  Flutter  Speeds  and  Frequencies  for  Three 
Configurations 


TARLE  2.2  Effect  of  Drag  on  Agreement  of  Modal  Analysis  with  Ref.  1 


Generalized 

P = 

0 • A 

P = 

0 . 00A 

Ampl i tude 

Phase 

Amplitude 

Phase 

Displacement 

qw 

1 

qw 

2 

qw 

0 . 71276 

219.62° 

1.A515 

195. A6° 

0.00132 

21.95° 

0.6288 

217. 2A° 

0.0000271 

150.1 3° 

0.00AA9 

203. 35° 

J 

qw 

A 

0 . 00000A08 

1A7 . 28° 

0.00067A 

-1 0 . 9A  ° 

’“5 

0.00000213 

163. A7° 

0.000189 

157. Al° 

\ 

1.0 

0° 

1.0 

0° 

'^2 

0.01503 

222.02° 

0.01A5A 

1A7. 11° 

0.000808 

215.78° 

0.01786 

210. 39° 

\ 

0.000773 

221. 31° 

0.00930 

210.13° 

\ 

0.000186 

217.90° 

0.00A62 

211.73° 

TABLE  2.3  Flutter  Mode  Shapes  for  Low-  and  High-Aspect-Ratio  Examples. 
(M  = AO.  ta  = 0.25,  A = 0.1.  S = 0.1.  C = 0.04) 
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;i 


i 


n 

mm 

C = 0.04 

UD 

% Difference 
From  n = 5 

UD 

% Difference 
From  n = 5 

i 

4.58288 

2.044% 

3.90105 

1.999% 

2 

4.48660 

-0.0999% 

3.81863 

-0.156% 

3 

4.49174 

0.0144% 

3.82465 

0.0017% 

4 

4.49067 

-0.0094% 

3.82403 

-0.0144% 

5 

4.49109 

— 

2.32458 

— 

(Ref.  1) 

4.66 

3.63% 

3.96 

3.42% 

TABLE  2.4  Comparison  of  Ref.  1 Results  With  (2-29)  Solutions  for 
Divergence  Speeds  in  the  Presence  of  Drag.  (P  = 0.004, 
M = 40.,  i = 0.25,  A = 0.1).  (For  C = 0,  Up  = 11.1072) 


r' 
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FIGURE  2-1 


FIGURE  2-2 
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Physical  Origin  of  Drag  Coupling  in  Eq.  (2-la) 


Physical  Origin  of  Drag  Coupling  in  Eq.  (2-lb) 
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Chapter  III 

DEVELOPMENT  OF  A GENERAL  SYSTEM  OF  EQUATIONS  FOR 
VERTICAL  BENDING/CHORDWISE  BENDING/TORSION  (VCT) 

A.  Introduction 

Chapter  II  includes  drag  loads  in  transverse  bending  and  torsion 
but  allows  no  chordwise  deflections.  In  other  words,  it  assumes 
infinitely  large  bending  rigidity  in  the  fore-and-aft  direction.  In 
this  chapter  we  consider  chordwise  bending,  which  is  a more  complete 
model  of  the  true  physical  situation. 

A review  of  the  literature  for  work  concerning  VCT  motion  of 
slender  cantilever  beams  led  to  the  field  of  hingeless  helicopter  rotor 
stability  analysis.  The  structural  modeling  of  a hingeless  rotor  is 
essentially  the  same  as  desired  here,  except  that  the  wing  has  no  rota- 
tional velocity.  Accordingly,  if  an  adequate  model  for  a hingeless 
rotor  can  be  found,  it  can  be  adapted  for  the  cantilever  wing  by  removing 
the  inertial  and  centrifugal  tension  effects  arising  from  rotation. 

Much  of  the  work  pursued  in  the  hingeless  rotor  field  includes 
simplifications  which  either  eliminate  or  restrict  one  of  the  three 
types  of  deformation,  often  torsion,  in  an  effort  to  reduce  complexity. 

This  leaves  a relatively  small  body  of  work  that  treats  the  full  elastic 

« 

problem.  A well-known  system  of  linear  partial  differential  equations 
for  coupled  elastic  torsion  and  bending  of  twisted  nonuniform  rotor 
blades  is  that  developed  by  Houbolt  and  Brooks  (Ref.  6);  the  initial 
effort  to  develop  a VCT  system  for  the  cantilever  wing  centered  on 
adapting  this  formulation.  In  the  course  of  this  work,  however,  the 
elastic  coupling  terms  were  found  to  be  insufficient  to  account  for  the 
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drag  loads  in  transverse  bending  and  torsion  embodied  by  the  terms 


containing  in  (2-1).  It  then  became  apparent  that  this  drag 

coupling  is  actually  a nonlinear  structural  bending-torsion  effect  and 
that  the  new  system  of  equations  to  be  developed  should,  of  necessity, 
retain  all  nonlinear  elastic  coupling  terms  having  the  same  order  of 
importance  as  these  drag  coupling  terms. 

Further  search  led  to  the  system  of  nonlinear  equations  for  twisted 
nonuniform  rotor  blades  derived  by  Hodges  and  Dowell  (Ref.  7).  This 
work  involves  development  of  a more  complete  strain-displacement  rela- 
tion than  that  of  Ref.  6,  which  is  necessary  to  obtain  the  elastic 
bending-torsion  coupling  terms  that  produce  the  desired  drag  coupling 
effect.  The  equations  in  Ref.  7 are  valid  for  straight,  slender, 
homogeneous,  isotropic  beams  undergoing  moderate  displacements,  accurate 
to  second  order  in  the  sense  of  a restriction  that  squares  of  bending 
slopes,  twist,  and  airfoil  chord  and  thickness  divided  by  wing  semispan 
are  small  with  respect  to  unity. 

Although  the  final  form  of  the  equations  presented  in  Ref.  7 might 
appear  to  be  immediately  adaptable  to  the  present  case  by  setting  the 
rotor  rotation  frequency  to  zero  and  removing  the  effects  of  pretwist, 
this  is  not  entirely  true.  One  important  assumption  by  Hodges  and  Dowell 
required  that  the  ratio  between  the  transverse  and  chordwise  beam  bending 
stiffnesses  be  a quantity  of  order  one.  Whereas  this  is  a standard 
feature  of  helicopter  blade  construction,  it  usually  does  not  hold  for 
conventional  aircraft  wings  of  any  aspect  ratio.  As  a result,  the  Ref.  7 
derivation  has  been  carefully  retraced  for  the  specific  case  of  a non- 
rotating cantilever  wing  having  arbitrary  bending  stiffness  ratio.  The 
development  is  outlined  in  the  following  section. 


i! 


1 
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B.  Development  of  the  Nonlinear  Equations  of  Motion 


Derivation  of  the  VCT  model  for  the  uniform  cantilever  wing  is 
presented  here  in  abbreviated  form.  Important  modifications  of  the 
Ref.  7 derivation  are  emphasized,  but  duplicate  manipulations  are  only 
briefly  described.  The  notation  and  coordinate  system  for  the  canti- 
lever wing  is  used  exclusively. 

The  basic  ordering  scheme  presented  in  equation  (4)  of  Ref.  7 is 
retained.  One  exception  is  that  spanwise  warping  of  the  cross  section 
due  to  twisting,  represented  by  a warp  displacement  function  which  is 
a solution  of  the  Laplace  equation  over  the  cross  section,  is  entirely 
neglected  here.  This  assumption  is  made  on  the  premise  that  a typical 
aircraft  wing  airfoil  section  would  have  a sufficiently  small  thickness 
that  warping  effects  would  be  negligible  within  the  second-order  frame- 
work. The  ordering  scheme  is  applied  to  the  energy  expressions  encoun- 
tered in  the  variational  derivation  of  the  equations,  to  determine  which 
terms  should  be  retained  and  which  discarded. 

The  nonlinear  strain  displacement  relations  developed  in  Ref.  7 
have  been  carefully  examined  in  the  context  of  the  present  problem,  and 
they  are  found  to  apply  without  modifications.  These  relations  are 
derived  from  an  exact  transformation  between  the  deformed  and  undeformed 
coordinate  sytems  and  originate  from  the  classical  definition  of  strain 
of  Novozhilov  (Ref.  19)  which  is  based  upon  increments  in  the  deformed 
coordinates.  After  approximation  to  second  order  consistent  with  all 
assumptions,  the  final  form  in  terms  of  engineering  strain  and  in  the 
present  notation  is 
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k 


A 


f 

,2  ,2 

u u 

£vv  = u'  + -y-  + -y u"  Kcos<t>  + nsin<t>] 

yy  y z z x 

- u"[-  Csin$  + r|cos<(>] 
z 


(3-1)  { 


= n4>  * 


These  can  be  recognized  as  identical  to  equations  (24),  (25),  (26),  and 

(27)  of  Ref.  7 after  the  warp  function  and  pretwist  angle  have  been 

eliminated.  The  displacements  u , u , u of  the  elastic  axis  and  the 

x y z 

principal  coordinates  £ , n of  the  cross  section  are  illustrated  in 
Figure  3.1. 

Development  of  the  equations  using  the  indirect  method  of  the 
calculus  of  variations  is  based  upon  Hamilton's  principle,  which  may 
be  stated  in  the  form 


(3-2)  / [ S(U-T)  - 6W]dt 

1 


The  equations  are  obtained  by  combining  expressions  for  the  first 
variation  of  strain  energy  6 U , kinetic  energy  ST  , and  virtual  work 
of  external  forces  Si'J  . 

The  first  variation  of  strain  energy  in  appropriate  form  for  the 
standard  Euler-Bernoulli  beam  uniaxial  stress-strain  relationship  is 


(3-3)  6 W = / JT  (a  6e  + a _ Se  c + o Se  )d£  dn  dy 

o area  yy  yy  yC  y$  yn  yn  s 


The  first  variation  of  the  engineering  strains,  expressed  in  terms  of 
displacements,  is 
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6u'  + u'  Su'  + u'  5u' 
y xx  z z 

- [f,c.os<f>  + rising ] (<Su"  - u"  Si^) 

X z 

- [£sin<J>'+  ncosij)]  (6u"  + u"  ) 

z x 

n 6 <j)' 

< 6 4>' 


The  stresses  are  simplv 


(3-5) 


0 

= Er 

yy 

yy 

0 _ 

= Gf  »- 

vf. 

yt 

0 

+ Ge 

y r> 

yn 

Substituting  (3-1),  (3-4),  (3-5)  into  (3-3)  yields 


( 3-h) 


$11  = JT  (E fu'_  + -y-  + — - u"({;cos<t>  + rising) 

area 

- u"(-£sln(|)  + ncoscjO  1 [(Su'  + u'cSu'  + u ’ (Su' 

Z it  V X X z z 

- (Ccosifi  + rising)  ((Su"  - u"  (Si^) 

X z 

- (-f,sin<f>  + ncoscb)  (iSu"  + u"  (S<J>)  ] 

Z X 

+ G(n?  4'  (S<K  + £2  4**  (S(f> * 1 kin  dr.  dy 


Rearrangement  of  (3-h)  bv  grouping  of  terms  having  the  same  virtual 
displacements  leads  to 
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(3-7) 


SU  = fQ  {Vy(6u-  + u^  6ux  + u’  6u')  + Sy  <5*' 
+ [-M  cos*  + M sin*  ] (<5u"  - u"  6*) 

£ X X 7. 


+ [Mz  sin*  + Mx  cos*](<5uj;  + 6*)}dy 


The  stress  resultants  and  moments  formed  In  (3-7),  which  act  on 
the  deformed  wing  as  illustrated  in  Fig.  3.2,  are  defined  as  follows 


(3-8) 


VV=  * °yy 

area  77 

,2  ,2 

u u 

= EAw { uy  + + ~ - eA(u^  cos*  - u^’  sin*)} 


where  A^eA  = JT  Z dn 


and  ff  n d£  dn  - 0 , by  definition 


M_  - H S or  dC  dn 
z vv 


(3-9) 


.2  ,2 

u u 


= El  (-u"  cos  + u"  sin  ) + EA  e.(u'  + ~ + ~) 
z x z w A v 2 2 ’ 


where  H ff  £2  d£  dn 


and  ff  £ n d£  dn  = 0 


\ = ~ n °vv  d^  dr^ 
area  77 


(3-10) 


EIx(ux  sin*  + u'z'  cos*) 


where  lx  = ff  n2  d£  dn 
area 


36 


(3-11) 


t' 


Sy  = SJ  (-£  a + n °xr)dC  dn 
area 

= GId  *' 

where  I,  = JJ  (£2  + n2)d£  dn 

d 

area 


At  this  point  in  the  development  the  moments  of  inertia  I , I 

are  first  introduced.  Hereafter,  it  must  be  recognized  that  I mav 

z 

assume  values  much  larger  than  I . The  terms  containing  the  quantity 
e^  , which  measures  the  offset  of  the  tensile  axis  (area  centroid)  from 
the  elastic  axis,  will  be  dropped.  The  basis  for  this  simplification 

is  that,  in  the  final  modal  formulation  of  the  equations,  the  retention 

EA  e . 2 
w A 

of  e^  terms  ultimately  leads  to  a dimensionless  parameter,  — — — , 

x 

which  appears  only  as  a small  quantity  added  or  subtracted  with  unity. 

Since  it  will  not  significantly  influence  dynamic  stability,  hereafter 

e.  = 0 will  be  assumed. 

A 

After  appropriate  integrations  by  parts  within  equation  (3-7),  the 
final  form  for  the  first  variation  of  strain  energy,  including  boundary 
terms,  is 


(5  U = f {[-(V  n ' ) * + (-M  cos<J>  + M sin<}')"]<$u 
o y x z x x 

- (V  ) ' 6u  + f—  ( V u')'  + (M  sin<t>  + M cosi}>) " l<$u 
y y y z z x z 


(3-12) 


+ [-(S  )'  - u"(-M  cos<J)  + M sin<J>)  + u"(M  sin<J>  + M cosi}')  ]<S<J>}dy 

V Z Z X X Z X 


+ V 6u  + uf  V 6u  I + u 1 V 6u  I ^ + S 

y vo  xy  x o z v z 1 o v 1 


o 


p p 

+ [-M  coscj'  + M sin(j>]6u'  | - (-M  cosi  + M sin<t>]'5u  | 

Z X XO  Z X xo 

+ [M^  sin<(i  + M^cosif  ]6i/ | ^ - [M^sin^  + M^cos^] '6u^  | ^ 
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Formulation  of  the  first  variation  of  the  kinetic  energy  is  greatly 
simplified  relative  to  the  Ref.  7 derivation  because  there  are  no  iner- 
tial effects  introduced  by  rotation  of  the  helicopter  blade.  Since  the 
procedure  is  straightforward  and  well  described  in  Ref.  7,  the  details 
of  forming  the  kinetic  energy  in  terms  of  displacement  velocities, 
taking  the  first  variation,  integrating  by  parts  over  time,  and  expressing 

the  resulting  form  of  <ST  in  terms  of  time  derivatives  of  u , u , 

x y 

u^  , and  4 are  omitted  here.  After  the  ordering  scheme  has  been 
applied,  the  form  of  the  first  variation  of  kinetic  energy,  with  terms 
retained  to  second  order  and  corresponding  to  equation  (52)  of  Ref.  7,  is 

&T  - JT  p (-  [ii  + (|>(-£sin4  + ncos4)]6u 
0 s x x 

area 


- [ii  - <f>(£cos4  + nsin<f>)  ]6uz 

(3-13) 


- [itx(-£sin<J>  + Ccostfi)  - uz(£cos<J)  + nsin<ji)  ] 6<Ji 

9 2 

- $[(-f,sin4  + 00034))“  + (£cos4  + nsin4)  l<$4)dn  d£  dy 

As  noted  in  Ref.  7,  the  last  term  in  this  expression  is  by  defini- 
tion a third  order  term,  but  it  is  retained  in  order  to  include  torsional 
inertia  in  the  torsion  equation. 

The  final  form  of  the  first  variation  of  kinetic  energy  is  obtained 
by  integration  over  the  sectional  area,  which  leads  to  the  following 
definitions: 


1 
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(3-14)  / 


m = If  p d?  dn 

s 

area 


a = If  P l d£  dn 
e s 


area 


If  Pg(£  + n )d£  dn 


area 


with  If  p n d^  dn  = 0 
area 


The  final  form  is 


(3-15) 


6T  ~ I {(-mu  + s </isin</>)  6u  + (-mu  + s $cos<J>)6u 

O X G X Z G Z 


+ [s  ii  sln</>  + s ii  cosP  - J<6]6<{>}dy 
ex  ez 


The  virtual  work  of  the  applied  loads  is 


(3-16)  60/  = / (F  6u  + F 6u  + m 6$)dv 

o x x z z y 


Clearly  the  drag  loads  will  now  enter  the  equations  of  motion  in  the 

same  manner  as  the  lifting  airloads,  in  contrast  with  (2-1),  since 

F and  F will  consist  of  lift  and  drag  force  components. 

X z 

Application  of  Hamilton's  principle  using  (3-12),  (3-15),  and 
(3-16),  together  with  (3-8),  (3-9),  (3-10),  and  (3-11),  results  in  the 
following  quantities  being  required  to  vanish: 


6u  terms 

— y 


(3-17a) 


= T'  = 0 
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<$u  terms 
— x 

{-El  (-u"cos(t>  + u"sin<j))cos<t> 
z x z 

{ 3- 1 7b ) + El  (u"sin4>  + u"cos4>)sin<i>}"  + mu 

x x z X 

— s 4>  sin$  = 0 
e 


6u  terms 
— z 


{El  (-u"  eos6  + u"  sin4>)sin<|) 
Z X z 


( 3- 1 7 c ) + El  (u^  sin<{>  + u^'  cos<f>)cos<|>}"  + muz 


- s $ eosd)  - F =0 
e z 


6d>  terms 


(GI,  4>  * ) * - u"  [ — El  (-  u"  cos<J>  + u"  sin<{>)cos<{> 
d z z x z 


( 3- 1 7 d ) 


+ El  (u"  sin<{)  + u"  coz<t>)sinc{>] 
XX  z 


+ u"  f El  (-u"  cos6  + u"  sind)) sir»(j)  + El  (u"  sin<|>  +■  u"  cosi)>)cos())] 

•• 1 XX  7 


X z X 


- s u sin<t>  - s u cos<£  + J<f)  - M =0 
ex  e z y 


In  (3- 17a) , since  it  is  known  that  for  the  nonrotating  contilever 
wing  the  spanwise  tension  T will  be  zero  everywhere,  the  expression 
can  be  integrated  leaving 
u’2  u'2 

(3-18)  u'  + -y-  + -j-  = 0 
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(true  only  with  = 0)  . This  constraint  indicates  that  the  only 
manner  by  which  axial  deflections  enter  the  problem  is  by  the 

purely  geometric  dependence  on  u^  and  u^  due  to  shortening  induced 
by  the  lateral  deflections. 

In  replacing  sin<{>  and  cosif  by  small  angle  approximations,  it 

would  be  consistent  with  the  ordering  scheme  to  approximate  the  cosine 

by  unity.  Due  to  the  possibility  of  a large  ratio  of  bending  stiffnesses, 

which  can  result  in  u^  being  large  compared  with  ux  , it  was  found 

necessary  to  keep  the  second-order  approximation  of  the  cosine  at  this 

step.  This  is  done  to  derive  properly  certain  elastic  coupling  terms 

involving  u and  u , while  maintaining  symmetry  in  the  inertial 

z x • 

and  stiffness  matrices  of  the  final  matrix  equations. 

Substitution  of 

COS<{)  ~ 1 - y- 

(3-19) 

sin<J>  ~ <p  4 

into  (3-17)  leads  to 

(3-20a)  [-  El  (-  u"(l-<t>2)  + u"(J>)  + El  (u"c}>2  + u"<J>)]" 

Z X Z XX 

+ mil  - s <6<b  = F 
x e x 

( 3-20b)  [El  (-  u"  + u"<J.)<|»  + El  (u”d>  + u"  (1-<J)2 ) ) 1" 

Z X Z XX  z 

+ mil'  - s $ = F 
z e z 
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(3-20c) 


(GI  <j>')'  - u"[-  El  (-u"(l-<J>2)  + u"<}))  + El  (u"4>z  + 

Q Z Z X Z XX  Z 


+ u"  [El  (-u"<f>  + u"<j>2)  + El  (u"<J>  + u"(l-((>2))] 

X Z X Z XX  Z Y 


- s u <t>  - s 11  +J$  = m 

e x e z T y 


Reorganization  of  terms  and  introduction  of  the  displacements  w 
and  v for  and  , respectively,  gives 


(3-21a) 


(3-21b) 


{'[EIz4>2  + EIx(l-cJ>2)  ]w"  - (EIz  - EIx)v"(()}" 


+ mw-s$  = F 
e z 


{-  (EI  - El  )w"cf>  + [El  (l-(t>2 ) + (f)2EI  ]v">" 
Z X z x 


+ mv-s(f)d)  = F 
e x 


GI.(J>"  - (EI  - EI  )[(w"2  - v"2)<t>  - v"w"  (l-2<p2 ) ] 
d z x 


(3-21c) 


sv(})  + sw  - J$  + m = 0 


The  underlined  terms  are  reasoned  to  represent  higher-order  effects  and 
are  dropped.  In  the  case  of  the  inertial  terms  coupling  chordwise 
bending  and  torsion  through  the  offset  , this  is  a third-order 

effect  relative  to  terms  like  mv  . The  final  form  of  the  three-DOF 
system  of  equations,  with  some  regrouping,  is 


( 3-22a) 


{eW'  - (EIz  - EIx)(v"<f>  - wV)}" 


+ mw  - s (6  = F 
e z 
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(3-22b) 


(El  v"  - (El  - El  )(w"<p  + y"<t>2)}"  + mv  = F 

Z Z X X 

( 3-22c)  GId4>"  - (EIz  - EI^)  [(w"2  - v£)<f>  - v'V] 

+ sw  = J$+m  = 0 
e y 

This  system  is  elastically  and  inertially  self  adjoint,  which 
assures  that  the  stiffness  and  inertia  matrices  eventually  assembled 
during  modal  analysis  will  be  symmetric.  The  terms  containing 
(EI^  - EI^)  represent  the  nonlinear  coupling  between  the  torsion  and 
bending  degrees  of  freedom;  all  remaining  stiffness,  inertial,  and 
applied-load  terms  represent  the  same  familiar  forms  encountered  in 
linear  beam  theory.  Terms  in  (3-22)  which  are  underlined  do  not  appear 
in  the  Ref.  7 equations  and  are  retained  here  as  a result  of  the  absence 
of  a restriction  on  the  bending  stiffness  ratio  EI^/EI^  . Strictly 
speaking,  when  this  ratio  is  large  compared  with  unity,  the  single 
underlined  terms  will  increase  in  relative  importance  while  the  double 
underlined  terms  are  negligible  in  magnitude. 

The  nonlinear  equations  of  motion  (3-22)  are  next  adapted  to 
permit  analysis  of  stability  about  a steady-state  deflected  position 
due  to  an  equilibrium  lifting  flight  condition,  which  could  be  level 
flight  or  a steady  pullout  at  a high  load  factor.  Small  time-dependent 
perturbations  about  the  equilibrium  operating  condition  are  introduced 
by  expressing  the  deflections  w,  v,  and  <J>  in  terms  of  steady-state 
equilibrium  deflections  w0  ’ v0  ’ and  <J>0  anc*  small  unsteady  perturba- 
tion quantities  w^  , v^  , and  <(>  : 


A3 


1 


(3-23) 


' w(t)  = w + W_  (t) 
o 1 

v(t)  = + v1(t) 

t 4>(t)  = <J>o  + <t>1(t) 


First,  the  steady  equilibrium  deflections  only  are  substituted 
into  (3-22)  to  obtain  a nonlinear  system  of  equations  for  the  equilib- 
rium solution.  The  resulting  nonlinear  steady  system  is 


( 3-24a) 

{El  w” 

- (El 

- El  )(v"4  - w"<J)2)}"  = F 

X o 

z 

X o o o o z 

o 

( 3-24b) 

{El  v" 

- (El 

- El  )(w"4>  + v’V)}"  = F 

z o 

z 

X O 0 O O X 

o 

(3-24c) 

GI  ,<p"  - 
do 

(EIz  - 

El  ) [(w"2  - v"2)<t>  - v"  w"]  + m 

X o o o O O V 

Appropriate  steady  lifting  aerodynamic  loads  are  inserted  for  L and 

w 

o 

M.  , and  the  assumed  steady  drag,  which  entered  equations  (2-1)  in  an 
vo 

entirely  different  manner,  is  introduced  through  L 

o 

Next  a linearized  system  of  equations  in  the  time  dependent  small 
perturbation  deflections  is  obtained  by  substituting  (3-23)  into  (3-22), 
subtracting  the  nonlinear  equilibrium  equations,  and  discarding  products 
of  the  perturbation  quantities.  The  linearized  unsteady  equations  of 
motion  are 


(3-25a)  {El  w','  - (El  - El  ) [<t>  vV  + *.v"  - <bzv"  - 2<J>  w"<t>.  ]}" 
x 1 z x o 1 1 o o 1 o o 1 


+ mw,  - s <5,  - F =0 
1 el  z 1 
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r 


(3-2 5b) 


(El  v'.' 
z 1 


(El  - El  ) [4>  V + W"4 
Z x o 1 o 


<f>2Vl  + 2<f>  v"*.  ]}” 
o 1 o ovlJJ 


+ mv,  - F 


( 3-25c) 


GI -4>’’  - (El  - El  ) { [ 2w"w"  - 
d 1 z x 1 o 1 


2v' V'  ] $ + 

o 1 o 


fWo2  - Vo2l*l 


o 1 


v"w"} 
o 1 


+ s w,  - J<S  + m 
el  ^1  y 


= 0 


Analysis  of  stability  with  this  system  can  be  done  with  the  standard 
techniques  for  linear  systems,  but  first  the  coefficients  must  be  found 
for  a given  equilibrium  flight  condition  by  solving  (3-24).  The  loads 
appearing  in  this  equation,  including  the  chordwise  forces,  must  be 
expressed  as  linear  functions  of  perturbation  displacements  w^  and 
<)>j  . No  dependence  of  aerodynamic  loads  on  the  fore-and-aft  motion 
v^  will  be  considered  in  this  analysis. 


C.  Comparison  of  the  VT  and  VCT  Models 

The  VCT  model  (3-24)  and  (3-25),  although  different  in  appearance 
from  the  VT  model  (2-1),  in  fact  reduces  to  the  same  form  when 
T ■+■  00  and  the  steady  lift  is  zero.  To  illustrate  this,  first  imagine 
chat  the  steady  loads  applied  to  (3-24)  are 


(3-26) 


F =0 
z 

o 


< m =0 

' y 

J o 


F = D 

^ xo 


* 


j 

! 
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where  D is  the  same  assumed  drag  force,  constant  in  magnitude  and 
direction  along  the  span,  as  is  considered  in  (2-2).  The  solution  to 
(3-24)  immediately  gives  = <J>  = 0 , leaving  just  the  chordwise 

bending  equation  as 


(3-27)  El  v""  = D 
z o 


Substitution  of  xEI  for  El  in  (3-25)  together  with 

x z 

w = =0  gives 

o o 


( 3-28a)  El  w’j"  - (x-l)EI  (<t>v")"  + mw  - s <6  - F = 0 

x 1 x ^1  o 1 eYl  z^ 


(3-28b)  xEI  vV"  + mv  - F =0 
x 1 1 x^ 


( 3-28c)  GI.<J>V  - (t-I)EI  [-v"^4>-  - v'V.’l  + s w - J$,  + m 
d 1 xl  o 1 o 1 el  1 \ 


The  equation  for  v^  uncouples  and  can  be  disregarded.  Integration 
of  (3-27)  twice,  using  the  zero  shear  and  moment  boundary  conditions 
at  the  free  end,  gives 


(3-29) 


tei^v^  = d Js(£-y)‘ 


= - M 


where  (2-2)  has  been  used.  Substitution  of  (3-29)  into  equations 
(3-28a,  c)  yields 

(3- 30a)  El  wV"  + — (M  $,)"  + raw.  - s - F =0 
v x 1 x zTl  1 el  z. 
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f 


(3- 30b) 


x— 1 


GI  ,<t>"  - M (w"  + v"<p. ) + s w - J<6  + m +0 

d 1 x z 1 o 1 el  rl  y. 


Finally,  requiring  that  x -►  00  , with  the  consequent  vanishing  of  v , 

o 

causes  (3-30)  to  reduce  exactly  to  (2-1). 

For  the  more  general  case  of  a steady  lifting  condition  charac- 
terized by  w ^ 0 and  $ ft  0 , (3-24)  and  (3-25)  can  still  be  cast 
o o 

into  a form  similar  to  (2-1),  in  which  the  v^  displacements  become 
dependent  upon  w^  and  since  X -*■  00  . First,  the  nonlinear  steady 

equations  (3-24)  are  rewritten  with  xEI^  replacing  EI^  together 
with  some  rearrangement  of  terms  to  give 


( 3-31a) 

El  w""  - xEI  [v" 

x o x X o 

- w"(J>  )(J) 

O O ( 

( 3—  31b ) 

XEI  [v"  - (w"<f>  + 

x o x o o 

v'V)]"  ■ 

o o 

( 3—  31c ) 

GI  ,<J>"  - — XEI  [ (w"<(> 
do  T x 1 oTo 

- vM)wn 
o o 

o 

..2,  i 
t <p 
o o 


v 

• o 


Again,  as  X 00  the  elastic  bending  curvature  about  the  n principal 

axis  (Fig.  3.1)  of  the  airfoil  section  should  go  to  zero.  But  now,  the 

section  is  displaced  to  a position  fixed  bv  the  deflections  w , v , 

o o 

and  (J)^  of  its  elastic  axis,  and  the  true  elastic  curvature  about  the 

H axis  is  now  recognized  as  v"  - w"d> 

o o o 

In  equation  ( 3—  31b ) , which  is  the  fore-and-aft  equilibrium  equation, 

as  x becomes  large  and  v becomes  small,  the  term  v"<)>2  is  of 

o o o 

higher  order  and  can  be  neglected  relative  to  v"  . Removing  this  term 

o 

X- 1 

and  letting  — — -*•  1 leaves 
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(3-32) 


TEI  [v"  - w"4>  ]"  = F 
xl  o o o x 


Integrating  this  equation  twice,  using  the  shear  and  moment  boundary 

conditions  on  v and  w at  the  free  end,  leads  to 
o o 


(3-33) 


TEI  [v"  - w"<p  ] = F ‘i(Z-v)2  = - M 
X O O O X z 

O o 


Here  it  is  assumed  for  convenience  that  F is  constant  along  the 

x 

o 

span;  of  course  the  integrations  could  as  well  be  performed  for  any 

known  spanwise  variation  of  F . M is  similar  to  M of  (3-29), 

o o 

except  that  it  represents  the  bending-torsion  coupling  effect  upon  the 

moment-curvature  relation.  As  T °°  the  quantity  (v"  - w"d>  ) must 

o o o 

go  to  zero  according  to  (3-33).  In  the  limit  the  v^  deflection  becomes 
dependent  upon  wq  and  4>q  as  a result  of  bending-torsion  coupling. 

The  two  remaining  independent  equations  in  (3-31)  are  already 
arranged  so  that  the  coupled  curvature  quantity  in  (3-33)  can  be 
recognized.  Sutstitution  for  this  quantity  leads  to 


(3-34) 


EIw""  + — [M  <f>  ]"  = F 

O T Z T O Z 

O 'o 


GI  ,(f>"  - — [M  w"  - v"2<{>  1 + m =0 
do  T z o o o y 
o ' o 


As  T -*■  00  , with  v”^p  recognized  as  a higher-order  term,  the  nonlinear 
steady  equilibrium  equations  finally  become 


(3-35a)  EIw”"  + [M  <J>  ]"  - F = 0 
o z o z 

o o 
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(3- 35b)  GI  , 4>"  - M w"  + m = 0 
do  z o V 
o ■'o 


These  resemble  (2-1)  with  dynamic  terms  removed. 

A final  exercise  is  to  demonstrate  that  the  linearized  unsteady 

perturbation  equations  (3-25)  also  reduce  in  a systematic  fashion  to  a 

VT  system  when  T -*■  00  . When  (3-25)  are  rearranged  with  tEI^ 

replacing  EI^  and  the  combination  of  steady  deflections  in  (3-33) 

replaced  by  M , the  linearized  system  becomes 
'o 


(3- 36a) 


( 3-  36b ) 


Eixwr + [mz  V - tetx  v - <w  - wv 


+ mw,  - s - F =0 
1 e 1 zt 


tF.I  (vV  - — (4>  w"  + w"4>  + <b2v"  + 2<(>  v'VMl" 

x‘  1 T VHo  1 oHl  Yo  1 o oH  1 


+ mv , - F =0 
1 


GI.<j>V  + TEl  — [(vV  - , w"  - w"i)w"  + 2v”4>  v."  + V"2<K  ] 

d 1 x T 1 11  o 1 o o o 1 o 1 


- M w"  + s w,  - Jif,  + m =0 
z 1 el  ’1  v. 
o 1 


The  underlined  terms  are  of  higher  order  for  t > 1 and  are  discarded. 

T—  | 

When  T is  sufficiently  large  that  — — ~ 1 , equation  (3- 36b)  becomes 


(3-37)  tEI  [v"  - ^ w"  - w"<|)  1"  = F - mv 
x 1 o 1 o’  1 x(  1 


(3- 36c) 
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The  quantity  $ , defined  as 

(3-38)  8 = v"  - <|>  w"  - w"tj>  , 

1 o 1 o 1 

is  recognized  as  the  time-dependent  curvature,  in  the  linearized 
perturbation  deflections,  about  the  n axis  of  the  steady-state  deformed 
cross  section.  Thus  it  is  expected  that  as  t •*  09  the  curvature  8 
must  approach  zero. 

The  dynamic  linear  moment  M is  defined  as 

Z1 

(3-39)  M = - XEI  8 
Z1 

Equations  (3-37)  and  (3- 36a,  c)  can  be  restated,  after  use  of  (3-38) 
and  (3-39)  together  with  T -*■  00  , as 

( 3- 40a)  El  w',"'  + [M  <J>,  ]"  + [M  <Jt  ]"  + mw  - s $ - F =0 

x 1 z 1 z , o 1 el  z , 

o 1 1 

(3- 40b)  M"  = F - mv, 

Z1  X1  1 

( 3-40c)  GI  ,</>''  - M w','  - H w”  + s w - J$,  + m =0 
d 1 z 1 z,  o el  1 v, 

o 1 1 

Since  8 0 , the  acceleration  term  in  the  second  equation  could  be 

expressed  entirelv  in  terms  of  w and  accelerations  hv  working 
from  integration  of 

(3-41)  v'j  = (J^w Y - w^t 

In  this  Linearized  unsteady  system,  M (v,t)  represents  the  moment  at 

Z1 

station  y due  to  the  Instantaneous  chordwlse  inertial  loads  and  applied 
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linearized  unsteady  drag  forces  acting  outboard  of  the  station.  Because 

of  the  dependence  of  upon  w^  and  , equation  (3-40b)  cannot 

be  uncoupled  to  allow  for  a separate  calculation  of  M"  . Thus,  even 

Z1 

though  v^  can  be  eliminated,  (3-40)  still  involves  three  unknowns 

w , 0 , and  M 

11  z^ 

The  purpose  of  this  section  has  been  to  demonstrate  analytically 
the  connection  between  VCT  representation  of  the  cantilever  wing,  given 
by  (3-24)  and  (3-25),  and  the  VT  system  given  in  (2-1)  and  used  in 
Ref.  1.  It  can  be  concluded  that  the  VCT  model  remains  valid  for 
arbitrarily  large  bending  stiffness  ratio  T , and  that  the  drag 
coupling  effect  has  been  satisfactorallv  accounted  for  by  the  nonlinear 
elastic  bending-torsion  coupling  terms.  Since  the  forms  obtained  in 
(3-35)  and  (3-40)  are  but  special  cases  of  (3-24)  and  (3-25),  actual 
solutions  will  be  found,  using  the  latter  system  only,  for  practical 
values  of  t . 
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Chapter  IV 


1 


FLUTTER  VELOCITY  CONSIDERING  VERTICAL 
BENDING/CHORDWISE  BENDING/TORSION  (VCT) 

A.  Modal  Equations  for  Small  Oscillatory  Motions  About  a Steady-State 
Deformation 

The  model  for  torsion,  transverse  bending,  and  chordwise  bending 

of  the  uniform  cantilever  wing  developed  in  the  previous  chapter  is 

used  to  determine  flutter  velocitv  with  of  the  same  set  of  assumed 

' 

modes  used  in  (2-1),  Use  of  n assumed  modes  in  torsion,  n in  trans- 
verse bending,  and  n in  chordwise  bending  results  in  a set  of  3n 
modal  equations  in  terms  of  modal  generalized  displacements.  Since  the 
assumed  modes  satisfy  the  natural  as  well  as  geometrical  boundary 
conditions,  which  were  obtained  during  the  application  of  Hamilton's 
principle,  Galerkin's  method  can  be  employed  to  transform  the  equations 
into  algebraic  relations  in  the  generalized  displacements.  The  nonlinear 
steady-state  equilibrium  equations  (3-24)  become  nonlinear  algebraic 
equations  in  the  steady-state  generalized  displacements,  which  are 
solved  iteratively.  These  displacements  determine  the  coefficients  of 
the  linearized  unsteady  model  by  applying  Galerkin's  method  to  (3-25). 
Then  the  velocity  is  determined  for  which  simple  harmonic  motion  of 
this  system  is  possible  (neutral  dynamic  stability). 

The  steady  aerodynamic  loads  for  the  equilibrium  equations  (3-24) 
are  specified  in  terms  of  incompressible  strip  theory.  A typical  air- 
foil section  (Fig.  4-1)  has  its  zero-lift  line  inclined  to  the  free- 

! stream  velocitv  V bv  the  angle  a + <J>  (v)  . The  resultant  steadv 

o ' 

lift  L (acting  in  a direction  perpendicular  to  V ) and  the  moment 
a 

; o 

m^  are  given  by 
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A drag  force  D is  also  present.  These  forces  must  be  resolved  into 
the  axis  system  x , y,  z fixed  with  respect  to  the  wing  root.  The 
required  transformation  is 


(4-2) 


/ 

F = L cosa  + Dsina 
z a 
o o 

IF  = Dcosa  - L sina 
x a 


Assuming  sina  ~ a and  cosa  ~ 1 and  neglecting  the  z component  of 
drag  shows  that  the  steady  aerodynamic  forces  to  be  used  with  (3-23)  are 


F = 2ITpVzb  (a  + <p  ) 
o 


(4-3)  < Fx  = 2IIpV2b(C  - a2  - a<(>o) 


m = 2ITpV2b2A(a  + $ ) 

v y0 


The  drag  force  represented  hv  C is  constant  spanwise,  and  the 
definitions  in  (2-22)  have  been  adapted. 

The  steady-state  deflections  are  now  expressed  in  terms  of  the 
assumed  modes  (cf.  Appendix  A for  definitions) 


J 
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— 


(4-4) 


n 

z 

f 

w . 

< 

i=l 

i 

n 

z 

f 

q° 

i=l 

vi 

V 

n 

z 

i=l 

ql. 

Bending  mode  shape  functions  f and  f are  identical.  Substi- 

Wi  Vi 

tuting  (4-4)  into  (3-24)  and  (4-3)  and  applying  Galerkin's  procedure 
leads  to 


/E L./l(E  f""q°  )f  dy  - (El  -El  ){/[(E  f"  q°  ) (Z  f.  q!  )]"f  dy 
i=l  i i j U=1  U y V=1  YV  j 


(4-5-)  < - {VC  qw.)(V*  q*  )(E1f4»q;)1"fw.dy} 

1=1  i l V=1  V v U=1  y y J 

^ = 2IIpVzb[a/f  dy  + J1  { E f.  q°  )f  dy] 

ow.  o .-d). 

j i=l  Yi  ri  j 


.El  /A(E  f""q°  )f  dy  - (El  -El  ){/£[(Z  f q°  ) (E  f.  q°  ) ]"f  dy 
2 0 i=l  vi  vi  V z x ° i=l  wi  wi  V=1  V*v  vj 


(4-5b)  { + />Z  r q-  )(E  f qj  )(E  f q;  ) ] " f y dy> 
i=l  l i v=l  V u=l  u y 3 


n 

2IIpVzb[(C  - az)Tf  dy  - a /*(  E f . q.°  )f  dy] 

j 1=1  Yi  i j 


(1  < j < n) 


GI.  f ( Z f"  q°  )f.  dy 
/do  v.=1  4>i  s<£>.'  ^ 7 


- (El  -El  ){/  ( E f"  q°  )(  E f"  q ° ) ( E f , q°  ) f . dv 

2 x 0 i=!  wi  wi  V-1  y=i  V*y  *j 


(4-5c)  I - fh  Z f"  q°  ) ( Z f"  q°  )(  Z f q“  )f 
) 0 i=!  vi  vi  v=i  Vv  vv  v-1  *y  ^y  *j 


dv 


- JT(  Z fM  q°  )(  Z f"  q°  )f.  dy 
o v v w w a) 

° i=l  vi  vi  x>  =1  v v v j 


V = - 2IIpV2b2A[a  f*  dy  + £ ’(  l ^dy ] 
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The  indicated  differentiations  in  the  first  two  equations  are  carried 
out,  and  appropriate  integrations  by  parts  of  the  resulting  modal 
integrals  is  performed.  The  first  terms  of  each  equation  are  then 
expressed  in  terms  of  modal  natural  frequencies  as  in  (2-10)  and  (2-11) 
which  leads  to  the  form 


™ - -rr  (t-1)  E EH.  q°  q 

j j y=i  v=i  Y\>  v 


(4-6a) 


n n n 


+ ~n~r  (T-1)  E E E R q°  q°  q° 
y=l  v=l  i=l  J H yy  i 

n f 2B.a 

- 2nov!b«  z i - 2nPv2b«  - 

V=1  ’ 


(A-6b) 


Tm  u2  2.  q 


w.  2J 

J j 


(t-1)  E Eh  . q!  q° 

v-l  u-1  “W  V“p 


n n n 


7T  (T-1)  l £ E R.  q!  q!  q° 

* v-l  y=l  i=l  ^ v' 


v y i 


+ 2HpV2b£a  E 7.  q°  = 2JTpV2b2  - rpl  (C  - a2)] 
v=l  JV  _ Nj 


(l  < j £ n) 


,2  2 0 , X 


n n n 


(4-6c) 


J < i q°.  + (T-1)  E E E R q°  q°  q° 

<p.  2 ‘<J>.  2 3 ' . ....  yvii  Mw  Mw  ‘dr 

Tj  v=l  y=l  i=l  ^ J v y i 

El  n n n 

- ~y£  (T-1)  E E E R . . q°  q°  q? 

1 v-l  y-1  i-1  WJi  \ \ \ 

= 2IIpV2b22A  a 


B.  and  N.  are  properties  of  the  assumed  bending  modes  described  in 
Appendix  A,  and  I..  represent  the  modal  integrals  previously  encountered 
and  defined  in  (2-13).  The  nonlinear  bending-torsion  coupling  terms 
give  rise  to  two  new  forms  of  modal  integrals, 


56 


w 


(4-7)  v ■ so  u <;>c  «)f: 

j u v 

(4-8)  Rijyv  = fo 

Equations  (4-6)  can  be  arranged  in  the  final  format  used  for  computa- 
tions by  dividing  the  torsion  equation  by  ITpV2£b2  and  the  two  bending 
equations  by  IIpV2£b  , using  modal  natural  frequency  relations  in  (2-14) 
and  (2-15),  and  nondimensionalizing  with  the  parameters  given  in  (2-16) 
and  (2-26). 


(4-9a) 


MPi  qw.  (T-l)MPi 

nV  O.  _J_ a 

1 u2  b u2^ 


n n n 


n n 


Z EH  . q I ~ 
y=l  v=l  \ b 


y=l  v-l  i=l  iJyv^y  % b -I 


l l E R.  . q°  q°  ~ - 2 Z I.  q° 

i ili\)  K I *(j) 

V=1  J 

n n ^w 

E £ H q°  . 

, , vyj  v b 

v=l  y=l  J v 


TMPi  Mv. 


MPi 


Vs  t1  - <*-»  v 


y 


4B . 

J_ 

TIN. 

J 


(4-9b) 


n n n 

+ E Z Z R,M,n,  9*  9h 


v=l  y=l  i=l  ljpv  ^y  b 


4B. 


+ 2a  £ I,  q 

A \ \ ‘ 


v=l 


N,n 


(C-a2 


(i  < j < n) 


(4-9c)  / 


9 Mirv 

[hn2(j-h)  - A]q? 


n ° n ° 

MPi  ("  n n n qw  he 

+ (T— 1)  — ^ j I l l R . . — ~ —H  q° 

U [y-1  u- 1 1=1  WVJ1  b b \ 


n n n 
E E E R 
v=l  y=l  i=l 


„ o O 

9 9 


M _V  o y y 

ell  "MV ji  b b q(J>±  yb1  Hjvy  b b 


q q 

v he, 

V _M 


= 2A  a 


H(j-Ji) 
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A solution  to  (4-9)  can  be  found  for  a given  wing  configuration 
once  drag  C , root  angle  of  attack  a , and  dimensionless  speed  U 
have  been  specified.  The  solution  procedure  involves  Newton's  method 
and  is  described  in  Appendix  C. 

When  simple  harmonic  motion  is  assumed,  the  linearized  perturbation 
equations  of  motion  (3-25)  are  converted  into  linear  algebraic  equations 
via  the  same  steps  described  in  Chapter  II  in  deriving  modal  equations 
(2-18)  from  (2-1).  The  perturbation  displacements  are  expanded  in 
terms  of  n assumed  modes  and  generalized  displacements. 


w^V.t) 


= £ f (y)q  e 


itdt 


. , w . ‘w . 

i=l  l l 


(4-10) 


v-,  (y,t) 


£ f (y)q  e 
. , w , v . 

1=1  l l 


iwt 


(y,t) 


« Vy>\ e 


lU)t 


where  q and  q have  dimensions  of  length,  whereas  q,  are 
w.  v. 

ii  l 


dimensionless.  Anticipating  solution  for  flutter  boundaries  using  the 
V-g  method,  structural  damping  g is  included  by  introducing  a complex 
elastic  modulus,  and  the  simple  harmonic  airloads  are  assumed  to  have 
been  expressed  in  terms  of  unit  generalized  forces  as  in  (2-9).  With 
the  addition  of  the  fore-and-aft  bending  degree  of  freedom,  the 
generalized  forces  now  appear  as 


(4-lla) 

(Uj^n) 


fo  Fz1(<t’l’wl;t)fw1  dy 


-WbHljj  qM-^+  ye1"* 
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(4— 1 lb ) 


'o  Fx1(*rwi5t)fVj  dy 


Tw  n 
1 + E 0, 


= npuzb*4[  E Q . -t—  + l l). . ,,0  q,  Je 

i=1  J+n.i  b i=1  ,1+n , i+2n  ^ 


ia)t 


(1  < j < n) 

(4-llc)  /?  m ^ (^i  ,wi  * c>  ^(j,  dv 


n w . n 

= TIpu>2b4£[  E Q „ - -r1  + E Q „ q.  ]e1(Jt 

,1+2n, l b s1+2n,i+2n  J 


The  omission  of  columns  in  the  array  Q.  . having  n +1  i 2n 
reflects  the  assumption  that  there  is  no  dependence  of  unsteady  air- 
loads upon  fore-and-aft  motions  . 

After  Galerkin's  method  has  been  implemented  for  equations  (3-25) 
by  the  same  manipulations  required  to  develop  (2-17),  the  system  of 
3n  modal  equations  which  determine  linear  stability  about  the  steady 
equilibrium  deflections  appear  in  the  form 


/ {M  - n4  N4  MPi  z}— r^- 
jab 


n n n w 

+ E (Q.  - (x-l)MPi  Z E E R,,  q”  q°  } 

. , j, i a , iUiv  ‘<J)  b 

i=l  J p=l  v=l  u 


(4-12a)  J + UlT-imi  z £ V,  "J  >TT 

) i=l  M=1  •'  W 


q q 

n n v n n w 

+ E { (i-l)MPi  Z[  E H.  . ■ - 2 E E R q°  -j-H 

. , a , ip  1 b , , Vivii  ^ b 

i=l  U=1  U=1  V=1  4u 


V ~ SMlJi+Qjii+2ntV  =0 
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(4-12b) 


n n V 

Z <Q1+n  , + (T-l)MPi  Z E H q°  }-ri 
i-1  “ u-1  ulj  *y  b 


+ {M  - TJI^N^fPi  Z}  — 
3 a b 


n n n v . 

+ E { (i-l)MPi  Z E Z R q®  q°  }— • 
. i a , , ijyv  ^d>  H<p  b 

i=l  y=l  v=l  J ry 


a q 

n n w n n v 

+ E {(t-I)KPI  Z[  E H.  + 2 E E R q° 

, , a 1 . iM  b , Vjyi  4>  b 1 

i=l  y=l  J y=l  v=l  J y 


}tU  = 0 


(1  < j < n) 


<j-fn,i+2n 


q q 

n ”i  v n n 

E { (x-l)MPi  Z[  E H . . — r^-  - 2 E E R . . q°  —H] 

. , a , yij  b , , yivi  44>  b 1 

l=1  y=l  H ■’  y=l  v=l  • Vv 


w . 

s m r..  + q _ .}  ~ 

lj  j+2n, i b 


(4-12c 


q Q q 

j r\  n n n v v . 

) { + E { (x-l)MPl  Z[  E H — H-  + 2 E E R . ,q° 

\ . ' a , yij  b . , yivj4^  b J-*b 

: 1=1  y=l  J y=l  v=l  ■'  v 


„ o „ o 

q q 

n n v . n n 


n nn  v vnn  q q 

+ E {(-r-l)MPi  Z[  E ER  . J-r£  - E E R y wv, 

• 1 a , , , yvij  b b , _ yvij  — -r— 

i=l  y=l  V=1  M J y=l  v=l  J b b 


+ V2n,i«„N1  + ‘■’"V1  - ■ 0 


Here  and  are  t^le  raocial  integrals  defined  in  (4-7) 

and  (4-8),  which  appear  in  this  system  multiplied  by  steady-state 
generated  displacements.  The  cofficients  of  terms  which  couple  the 
steady  deflections  into  this  system  must  first  be  constructed  by  the 


indicated  summations  before  the  eigenvalues  can  be  determined. 


The  unit  generalized  forces  required  in  (4-12)  will  now  be  expressed 
in  terms  of  the  two-dimensional  incompressible  unsteady  aerodynamics 
adopted  for  (2-2)  and  (2-3).  It  will  be  necessary,  as  it  was  for  the 
steady  aerodynamics,  to  transfer  the  airloads  from  wind  oriented  axes 
to  the  coordinate  system  fixed  at  the  wing  root.  For  the  present  case, 
the  circulatory  part  of  the  unsteady  lift  will  be  assumed  perpendicular 
to  the  direction  of  the  freestream  velocity,  whereas  the  noncirculatorv 
portion  will  be  assumed  to  act  normal  to  the  chord  of  the  airfoil  sec- 
tion in  its  steady-state  deflected  position.  The  inclusion  of  unsteady 
leading-edge  suction  effects  will  be  considered  in  the  next  chapter. 

To  separate  the  circulatory  and  noncirculatorv  contributions  to 
the  lift,  L in  (2-2)  and  (2-3)  can  be  rewritten  in  the  form  (cf.  Eq . 
(4-126)  of  Ref.  5) 

w i(fi1 

(4-13)  La(w1,4>1;t)  = Hpb*oj2[-^-  + + a^e1^ 

w 2$' 

+ flpb  Vc(k)  [-  + -j-^  + ^ (%-a)<|)1]elut 


Here  the  first  term  on  the  right  represents  noncirculatorv  lift  I. 

aNC 

and  the  second  term  the  circulatory  portion  L . Referring  to 

ac 

Fig.  4-2,  these  terms  can  be  expressed  as  resultant  forces  in  the  xz 
frame  by 


(4-14) 


L cosa  + L cos<i> 
3C  aNC 


- L sina  + L sin<}> 
3C  aNC 
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Making  the  assumptions  cosa  ~ 1 and  costfi^  ~ 1 yields  the  identity 


(4-15)  Fz  = La(w1,<J)1;t) 


When  it  is  assumed  that  sina  ~ a and  sin<)>  ~ d)  , the  chordwise 

o o 


force  components  become 


(4-16)  F = npbVa  C(k)[^  ^ " pr  - ^ (*s-a)$1  ]eia)t 
+ npb3oo24>  [-;r  + ~r^  + ]elu)t 

ob  k i 


Preparatory  to  steps  that  lie  ahead,  it  can  be  seen  that  when  the 
unit  generalized  forces  are  determined  from  (4-11)  with  the  modal  series 
for  <J>o  , w^  , and  <J>  inserted,  modal  integrals  of  the  forms 


fo  fWl(?)fv  (y)f<t> (?)d? 

1 j 

fo  fv.(?)f4)  (?)d? 

i j 

will  be  encountered  as  a consequence  of  the  noncirculatory  contribution 

to  (4-16).  Since  these  integrals  do  not  occur  elsewhere  and  the  effect 

they  introduce  is  expected  to  be  minor,  the  noncirculatory  contribution 

to  F is  neglected  altogether.  This  simplification  has  the  result 
X1 

that  the  unsteady  noncirculatory  force  illustrated  in  Fig.  4-2  acts  in 

the  z direction  at  all  spanwise  stations. 

The  remaining  circulatory  contribution  to  F can  be  expressed, 

X1 

using  the  notation  in  (2-3)  and  after  some  manipulation,  as 
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(4-17) 


F = npbVa{(l-L  ) ~ + [L.  - h + £ + A(l-L  )]3>  e 

W D (J)  k W 1 


With  the  aerodynamic  loads  given  in  (4-17),  (4-15),  and  (2-2), 

where  m (<f>  ,w  ;t)  = m (w_  , 4>,  ; t ) , the  unit  generalized  forces  can 
l 1 3 11 

be  developed,  starting  from  (4-11),  by  means  of  the  same  steps  followed 


in  assembling  (2-22).  The  final  result  is 


f Q1,i 


Lw  (1=J> 

0 (ifl) 


j , i+2n  = - (L<|>  - AL  ) I, . 

w j 1 


.1+n,  i 


( 1-L  )ot  (i=j) 

w 

0 (i#j) 


(4-18) 


(11.1£n) 

(l<i<n) 


Q..  * [L.  - h + r-  + A(l-L  ,)]I.  . 

\l+n,i+2n  ip  k w ij 


Q, , = - (M  - AL  )I. , 

j+2n,i  w w ij 


Qi+2n,i+2n  - A(L^  + Mw)  + A2!.^  1 (i-j) 


0 ( 1 i ) 


■ 0 for  I < v < In 


■ ™ — 


B.  Solution  Procedure  for  Flutter  Velocity  of  a Lifting  Wins 

The  3n  modal  equations  (4-12),  with  (4-18)  inserted,  can  be 
expressed  in  matrix  form 

(4-19)  ( [Mg]  + [ Q ] ) {q  1 - Z[kJ{q}  = 0 

The  mass  and  stiffness  matrices  M and  K are  real  and  symmetric 

s s 

and  are  written  out  in  Fie.  4-3.  The  terms  in  K include  sums  of 

s 

products  of  modal  integrals  and  the  steady  equilibrium  generalized 

displacements  found  as  solutions  of  (4-9).  The  aerodynamic  matrix 

Q , whose  elements  appear  in  (4-18),  is  complex  and  nonsymmetric . 

Equation  (4-19),  therefore,  represents  a complex  eigenvalue  problem 

for  the  complex  frequency  parameter  Z . Its  solution  yields  damping, 

speed,  and  frequency  as  in  the  case  of  (2-25)  and  (2-27). 

The  logic  used  to  compute  neutral  stability  conditions  from 

(4-19)  is  diagrammed  in  Fig.  4-4.  The  primary  difficulty  encountered 

when  steady  deflections  are  introduced  is  that  a preliminary  estimate 

of  speed  must  be  made  before  the  eigenvalue  problem  can  be  solved. 

Steady  deflections  for  U are  used  to  generate  coefficients  in  K , 

e s 

and  unsteady  aerodynamic  loads  are  then  computed  for  large  enough 
reduced  frequency  that  the  eigenvalue  corresponding  to  the  aeroelastic 
mode  which  flutters  has  negative  (stable)  damping  g . Successively 
smaller  values  of  k ire  then  substituted  and  aerodynamic  terms 
•.  ■ a.  Pie  eigenvalues  ire  recalculated  until  positive  damping 

«nb- 


it 


is  repeated  with  a new  estimate  U until  the  speeds  U and  U_ 

e e F 

are  adequately  matched.  The  flutter  speed  determined  for  zero  lift 
serves  as  a good  first  estimate  for  . With  care,  close  agreement 

of  the  two  speeds  can  be  achieved  in  three  or  four  iterations. 

Three  assumed  modes  in  each  degree  of  freedom,  corresponding  to 
a system  order  of  9,  are  found  to  give  adequate  convergence  for  all 
cases.  The  model  integrals  were  computed  to  allow  n < 4 , a task 
which  required  numerical  integration  of  100  quantities  of  type  R ^ 
and  40  of  type  H. 


C.  The  Nonlinear  Elastic  Coupling  Terms 

The  need  to  retain  all  terms  in  (3-22),  including  third  degree 
nonlinear,  in  order  to  model  adequately  the  nonlinear  elastic  bending- 
torsion  coupling  mechanism  is  now  demonstrated  by  means  of  typical 
applications.  Of  course,  neglecting  higher-order  nonlinear  effects 
would  have  the  appeal  of  reducing  complexity.  For  example,  removal 
of  all  third-degree  nonlinear  terms  would  result  in  the  elimination 
of  all  terms  containing  quadruple  modal  integrals  R ^ from  both 
the  nonlinear  steady  equilibrium  system  and  the  linearized  dynamic 
stability  analysis.  In  order  to  examine  the  effects  of  such  approxi- 
mations, numerical  experiments  were  conducted  wherein  higher  degree 
terms  in  both  the  steady  and  unsteady  modal  equations  were  neglected. 

First,  static  deflections  are  considered.  Clearly,  if  the 
equations  do  not  adequately  represent  the  steadv-state  deflections 
of  the  winit  over  i real  let i<  ranee  of  1 if  tine  conditions,  then  anv 


approximation  cannot  be  expected  to  succeed.  The  nonlinear  steady 
equations  (3-23)  will  be  solved  at  three  levels  of  approximation: 


(1)  Linear  terms  only  retained. 

(2)  Second  degree  monlinear  terms  included. 

(3)  All  nonlinear  terms  included. 

The  linear  set  of  modal  equations  is  just  (4-9)  with  the  nonlinear 
terms  removed : 


- 2 
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Z I,  q! 

V=1  'V 


(4-20) 
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The  torsion  equations  are  now  seen  to  be  n uncoupled  and  immediately- 
solveable  relations  for  the  n generalized  displacements  q°  . The 


j 

results  permit  solution  for  the  bending  displacements,  thus  eliminating 
the  need  for  matrix  operations.  Since  all  elastic  coupling  terms  are 
absent,  deflections  w^  and  <|>  are  independent  of  the  stiffness 
ratio  T and  the  steady  drag  parameter  C . A very  important  conse- 
quence is  that  the  mechanism  by  which  drag  influences  divergence  is 
missing.  One  concludes  that  (4-20)  is  a particularly  poor  representation 
for  steadv  deflections  in  the  presence  of  steady  drag. 

Retention  of  the  second-degree  nonlinear  effects  from  (4-9)  involves 
nliling  to  (4-20)  all  terms  containing  the  triple  modal  integrals 
M . The  system  be<'*s»*  till  I v coupled,  and  solutions  are  found  bv 


3 


the  same  iterative  scheme  described  in  Appendix  C for  solving  the  full 
nonlinear  system. 


1 


The  behavior  of  the  second  degree  nonlinear  solution  is  compared 
with  the  linear  deflection  from  (4-20)  in  Figs.  4-5  and  4-6  . Here 
it  is  shown  how  the  deflection  in  the  first  vertical  bending  mode  and 
the  first  torsion  mode  vary  with  increasing  speed  for  a fixed  angle  of 
attack  a = 0.01  radions  of  the  wing  root.  The  parameters  M , P , 
ia  , A , and  T for  these  examples  correspond  to  the  idealization  of 
a sailplane  wing  discussed  in  Chapter  VI,  and  sufficient  modal  con- 
vergence is  assured  by  using  n=3  . Numbers  on  the  ordinate  of  Fig.  4-5 
show  the  actual  twist  in  radians  of  the  wingtip  due  to  the  first  torsion 

mode;  the  vertical  deflection  in  semichords  of  the  wingtip  due  to  q° 

W1 

is  just  twice  the  value  read  from  the  ordinate  of  Fig.  4-6. 

The  important  influence  of  drag  on  the  steady  deflections  is 
evident.  But  there  is  poor  correspondence  between  the  second  degree 
nonlinear  deflections  and  the  divergence  speeds  indicated  in  Fig.  4-5  , 
which  are  solutions  of  the  linear  static  stability  determinant  (2-28).  i 

This  discrepancy  is  most  pronounced  when  C=0,  and  the  rapid  divergence 
of  the  second-degree  nonlinear  deflections  for  all  values  of  drag 
suggests  trouble  with  the  second-degree  approximation. 

Figures  4-7  through  4-10  show  solutions  of  the  complete  non- 
linear system  (4-9)  for  the  same  wing  configuration,  together  with  the 
linear  results.  Correlation  with  the  divergence  predictions  appears 
to  be  excellent,  and  the  sudden  blowing  up  of  deflections  characteristic 
of  the  second-degree  nonlinear  solutions  is  not  encountered.  Since  the 


with  the  linear  solution  as  an  initial  estimate,  nonlinear  solutions 
which  do  exist  for  C=0  above  the  classical  divergence  speed  cannot 
be  obtained. 


The  full  nonlinear  solutions  conform  with  nonlinear  behavior 
expected  by  intuition.  For  zero  drag,  the  vertical  bending  and  torsion 
deflections  should  fall  below  the  linear  solution  owing  to  the  effec- 
tive increase  in  stiffness  "seen"  by  each  degree  of  freedom  due  to 
deformation  in  the  other.  This  effect  is  observed.  The  slight  rear- 
ward chordwise  displacement  (Fig.  4-9)  for  C=0  when  w and  4> 

o o 

are  large  comes  from  bending-torsion  coupling,  and  the  negative  contri- 
bution of  the  second  torsion  mode  (Fig.  4-10)  reflects  a redistribution 
of  elastic  twist  toward  the  wingtip  where  the  curvature  due  to  w^  is 
less.  Drag  alters  the  deformed  state  radically  and  causes  large 
displacements  at  much  lower  speeds. 

A clue  to  the  reason  why  the  second-degree  terms  alone  are  inade- 
quate is  found  by  looking  at  the  sensitivity  of  solutions  to  changes 
in  the  bending  stiffness  ratio  T . The  second  degree  nonlinear 
solution  for  bending  deflections,  given  in  Fig.  4-11,  reveals  increas- 
ingly poor  behavior  as  T is  increased.  Conversely,  the  complete  non- 
linear solution  behaves  as  intuition  would  anticipate,  becoming 
insensitive  to  changing  T as  this  parameter  grows  toward  00  (Fig.  4-12). 

In  Chapter  III  the  behavior  for  of  the  steady  nonlinear 

equilibrium  equations  was  analytically  investigated.  In  the  discussion 
following  equation  (1-13)  it  was  pointed  out  that  the  quantity 

v"  - * w"  represents  the  actual  curvature  about  the  principal  axis 

o o a 

if  < hordwlsr  hemline  for  the  deflected  tfrfoil  section,  ind  that  this 


quantity  should  vanish  as  T ■+  00  . In  this  limit  any  chordwise 

deflection  is  a geometric  consequence  of  elastic  coupling  between  <)> 

and  ; terms  which  contain  the  product  t(v^  “ <J>ow^')  approach  a 

finite  value.  Now,  referring  to  equations  (3-24),  it  can  be  seen  that 

in  the  second  degree  approximation  this  product  is  retained  in  the 

second  equation  but  not  in  the  vertical  bending  equation  (where  w"d)2 

o o 

is  dropped  but  v^4>0  retained).  Likewise  in  the  torsion  equation, 
w^2<f>o  is  dropped  but  v^w^  is  kept.  The  result  is  that  terms 
containing  Tv^  remain  and  blow  up  in  the  limit  T -*•  00  . Accordingly, 
the  conclusion  is  reached  that,  for  structures  representative  of  air- 
craft wings  for  which  T is  reasonably  large,  the  third-degree  terms 
containing  the  product  w^(J>o  must  be  retained.  It  may  be  added  that, 
although  this  elastic  coupling  effect  was  neglected  in  the  helicopter 
blade  equations  developed  in  Ref.  7,  in  that  case  the  ordering  scheme 
required  x to  be  on  the  order  of  unity.  Another  point  is  that  the 
remaining  third-degree  terms  in  (3-24),  which  contain  the  product 
v^$o  , actually  can  be  neglected  when  T is  large. 

Since  it  has  been  determined  that  third  degree  terms  cannot  be 
excluded  from  the  steady  nonlinear  equilibrium  analysis,  it  follows 
that  they  must  also  be  kept  in  the  linearized  unsteady  perturbation 
equations  (3-25).  The  original  dynamic  stability  equations  used  early 
in  this  investigation,  on  the  other  hand,  consisted  of  a linearized 
unsteady  perturbation  system  based  on  the  second-degree  nonlinear 
approximat ion.  That  Is,  thov  are  the  result  of  removing  all  terms 
containing  the  modal  integrals  R(^  from  (4-12).  Irregular  behavior 


*f  the  tint  ter  speed  with  increasing 


I attack 


vered , 


as  was  an  extreme  sensitivity  to  large  T which  could  not  be  justified 
physically.  Furthermore,  linear  steady-state  deflections  were  used — 
a particularly  unsuitable  approximation  for  inclusion  of  steady  drag. 

Examples  of  these  original  flutter  calculations  appear  in  Figs. 

4-13  and  4-14.  The  former  gives  a comparison  of  the  early  results 
with  the  full  analysis  by  (4-12)  for  the  sailplane  example,  showing 
how  the  simpler  approximation  differs  significantly  in  flutter  speed 
even  in  the  presence  of  moderate  steady  deflections.  Figure  4-14, 
based  on  a different  wing  configuration,  gives  an  idea  of  the  difficulty 
encountered  for  large  T when  using  the  simpler  analysis.  (Note  that 
the  steady  deflections  are  quite  small  in  view  of  the  fact  that  this 
example  is  a large-aspect-ratio  wing).  All  of  these  stability  boundaries 
abruptly  terminate,  at  which  point  the  eigenvalue  solutions  began  to 
behave  erratically.  In  contrast,  flutter  solutions  of  (4-19)  can  be 
obtained  for  arbitrarily  large  steady  deformations — indeed  for  tip 
deflections  well  beyond  practical  limits  of  material  linearity. 


just  vertical  bending  and  torsion. 


The  VCT  stability  equations  were  analytically  reduced  to  the  VT 
form  (2-1)  in  the  last  chapter  by  specifying  a = 0 and  x -*■  00  . This 
agreement  is  reflected  in  numerical  results,  as  evidenced  by  Fig.  4-15. 
For  an  exceptionally  large  magnitude  of  drag,  C = 0.04  , chosen  to 
magnify  the  importance  of  drag  coupling  and  hence  the  steady  deflections 
in  the  stability  analysis,  the  two  methods  are  used  to  compute  flutter 
speeds  for  the  same  nonlifting  configuration.  The  VCT  analyses  were 
made  for  T ranging  from  1 to  10,000. 

As  is  always  true  for  EI^  = EI^  , the  drag  coupling  effect 
vanishes  at  T = 1 with  flutter  speed  unaffected  by  drag.  At  the 
other  extreme,  for  x = 10,000  the  computed  flutter  speeds  for 
C = 0.04  differ  by  a mere  0.052%.  The  Ref.  1 flutter  speeds  calculated 
for  this  same  configuration  are  also  shown  and  agree  quite  well. 

When  one  examines  the  VCT  dependence  of  flutter  speed  on  x 
(Fig.  4-15),  the  effect  of  drag  is  apparently  insensitive  to  X over 
the  range  of  this  parameter  representative  of  practical  aircraft 
applications.  Hence  the  VT  system  (2-1)  does  not  suffer  by  its 
inherent  assumption  that  x -►  00  . 

Table  4.1  presents  a further  comparison  of  results  from  the  two 
approaches  for  a different  configuration,  in  which  modal  convergence 
is  emphasized.  Of  course,  since  X = 50  exact  agreement  cannot  be 
expected . 

Convergence  of  flutter  speed*  and  frequencies  and  also  mode  shapes 
for  the  VCT  solution*  Is  doruMS  t ed  In  Table  * . for  rero  '*f»-*dv 


lift  and  in  Table  4.3  for  a lifting  condition  specified  by  a = .01 
rad.  In  the  latter  case,  the  apparent  slower  convergence  is  related 
to  the  steady  lifting  deflections  (also  tabulated),  which  also  vary 
with  n . 


A 


UF,  VT 

UF,  VCT 

4.260823 

4.258457 

4.260879 

4.258351 

4.260882 

4.258335 

4.260889 

4.258336 

TABLE  4.1  Comparison  of  Modal  Flutter  Speeds  Computed  with  the 
VT  and  VCT  Systems.  (M  = 40 . , P = 0 . 4 , i =0.25, 
A = 0.1,  S = 0.1,  T = 50,  C = 0.04)  01 


n = 2 

n = 

3 

n = l 

U_.  = 3.3495 

F 

UF  = 

3.5329 

U-,  = 3.5315 

F 

= 0.68200 

“p" 

0.67468 

f2p  = 0.67276 

o 

qw 

1 

1.19055 

1.18142 

1.18074 

qw 

2 

0.01336 

0.01333 

0.01333 

<c 

- 

0.000817 

0.000818 

3 

qw 

4 

“ 

- 

0.000151 

O 

qv 

0.002270 

0.002296 

0.002305 

V1 

-0.000279 

-0.000284 

-0.000284 

<3 

-0.000033 

-0.000033 

-0.000005 

o 

\ 

0.0087100 

0.0087591 

0.0087656 

q‘2 

0.0000373 

0.0000915 

0.0000983 

0.0000017 

0.0000090 

“ 

- 

0.0000013 

AMPL. 

PHASE0 

AMPL. 

"i : 

, PHASE0 

1 

AMPL. 

PHASE0 

qw 

1 

2.4492 

227.45 

2.4765 

| 227.68 

2.4860 

227.76 

% 

0.52487 

-26.65 

0.50948 

| -26.27 

0.51068 

-26.23 

2 

1 

%, 

0.01637 

j - 6.54 

0.01612 

- 6.57 

3 

1 

% 

4 

“ 

1 

0.00250 

- 1.47 

qv 

0.84595 

- 4.77 

0.85914 

i - 4.89 

0.86569 

- 4.98 

V1 

1 

qv 

V2 

0.03071 

179.03 

0.02979 

I 180.02 

0.02916 

180.20 

qv 

0.00923 

I 172.67 

0.00532 

172.76 

V3 

\ 

- 

1 

0.00122 

172.56 

\ 

1.0 

0 

1.0 

1 0 

1.0 

0 

^2 

0.07178 

-21.65 

0.07798 

-20.83 

0.07947 

-20.73 

qA 

♦3 

0. 027<*  3 

| - 4.83 

0.02925 

- 5.06 

V 

• 

- 

0.00998 

- ft  .f»0 

ri  j 


[Mg ] = M 


“S*^ll  'ScAn 


'Hn 


C~I  J 


J 


ABBREVIATIONS  OF  SUMMATIONS  APPEARING  IN  [K  ] : 

s 

n n 

R.  = 2 2 R.  . ql  ql 

ij99  L V = 1 <Pm  9v 


R ...  = X 2 R . . 
wi</>J  /i  =1  v = ! vl/lJ  Vn 


o 

o S/v 


H,  . . = 2 H . . q , 

91  J ,.=i  viJ  9v 


n q. 

H.  . = I H.  . - v 

iwj  _ i i*'j  b 


FIGURE  4- 3(a)  Mass  and  Stiffness  Matrices  Appearing  in  Equation 
(4-19) 
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FIGURE  4- 3(b) 


FIGURE  4-4  Flow  Chart  for  Solution  Procedure  Using  V-g  Method 
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FIGURE  4-5  Deflection  of  the  First  Torsion  Assumed  Mode  for  « = .01 
and  Increasing  Speed,  Comparing  Linear  and  Second-Degree 


Nonlinear  Solutions. 
A * 0. 1 , S * 0. 1,  T = 


(M  = 9.4,  P = 0.01  , t 
25,  and  n = 3) 


0.25, 


FIGURE  4-6  Amplitude  of  the  First  Vertical  Bending  Mode.  (Same  Conditions  as  in  Figure  4-5) 


0.  1.  2.  3. 

u 


FIGURE  4-7  Amplitude  of  the  First  Torsion  Mode  for  a = 0.01  and 
Increasing  Speed,  Comparing  Linear  and  Full  Nonlinear 
Solutions.  (M  = 9.4,  P = 0.01,  i^  = 0.25,  A = 0.1,  S 
T = 25,  and  n = 3) 
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NONLINEAR 


FIGURE  4-9  Amplitude  of  the  First  Chordwise  Bending  Mode.  (Same  Conditions  as  in  Fig.  4-7) 
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Figure  4-10  Amplitude  of  the  Second  Torsion  Mode.  (Same  conditions 
as  in  Fig.  4-7) 
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FIGURE  4-11  Amplitude  of  the  First  Vertical  Bending  Mode  vs.  T,  Using  Second-Degree  Nonlinear 
and  Linear  Approximations.  (Remaining  Parameters  Same  as  in  Fig.  4-5) 


nonlinear 


FIGURE  4-12  Amplitude  of  the  First  Vertical  Bending  Mode  vs.  x.  Using  Full  Nonlinear  and  Linear 
Approximations.  (Remaining  Parameters  Same  as  in  Fig.  4-7) 


FIGURE  4-13  Flutter  Speed  as  a Function  of  Steady  Wing-Tip  Displacement  for  Two  Different 
Approximations  in  the  Analysis.  (M  = 9.4,  P = 0.01,  i = 0.25,  A = 0.1,  S = 


Chapter  V 


THE  DETERMINATION  OF  AEROELASTIC  MODES  FOR  ARBITRARY  VELOCITY 

A.  Incompressible  Strip  Theory  Airloads  for  Arbitrary  Motion 

The  procedure  for  determining  flutter  velocity  developed  in 
Chapter  IV  has  several  drawbacks.  It  requires  a matched  point  analysis 
in  which  an  estimated  velocity  used  to  calculate  steady  deflections 

has  to  be  iteratively  matched  to  the  lowest  calculated  flutter  speed 
U for  the  proper  structural  damping.  Consequently  the  solution  for 

r 

a stability  boundary  over  a range  of  lifting  conditions  can  be  lengthy. 
Furthermore,  intermediate  computations  have  no  physical  significance 
and  are  of  limited  qualitative  value.  The  behavior  and  degree  of 
stability  of  individual  aeroelastic  modes,  which  becomes  more  inter- 
esting with  the  addition  of  the  fore-and-aft  bending  degree  of  freedom, 
has  proven  to  be  difficult  to  deduce  from  the  V-g  solutions.  The 
only  quantitative  information  available  pertains  to  the  neutral  stability 
conditions  found  for  the  mode  which  experiences  flutter. 

We  would  like  to  have  physically  meaningful  information  regarding 
dynamics  of  the  system  at  any  desired  speed.  That  is,  we  would  like 
to  know  the  complex  eigenvalues  of  the  aeroelastic  modes  at  subcritical 
and  supercritical  velocities.  Obviously  a major  drawback  of  the  V-g 
method  is  its  dependence  upon  simple  harmonic  air  loads. 

Several  solution  procedures  were  studied,  which  replace  the 
Fourier  transformation  with  respect  to  t by  Laplace  transformation 
with  respect  to  t and  obtain  at  least  an  approximation  to  the  modal 
stability  below  and  above  the  flutter  velocity.  The  p-k  method,  a 
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British  flutter  analysis  technique  (Ref.  8),  uses  the  same  simple 
harmonic  airloads  but  assumes  that  these  loads  are  approximately  com- 
plex eigenvalues  s = a + ioj  , where  a)  is  the  frequency  used  in 
calculating  the  airloads  and  |o|  <<  jio|  . 

Another  approach,  commonly  applied  in  helicopter  blade  stability 
analysis,  involves  use  of  quasisteadv  aerodynamic  theory  by  assuming  k 
is  small  enough  to  permit  C(k)  ~ 1 . Laplace  transformation  of  the 
system  then  yields  only  linear  and  quadratic  terms  in  s . Linear 
matrix  techniques  can  then  be  used  in  determining  the  roots.  This  level 
of  approximation  neglects  entirely  the  effect  of  the  unsteady  wake  upon 
the  circulatory  airloads  and  is  not  appropriate  for  the  magnitudes 
of  reduced  frequencies  observed  in  many  flutter  calculations  by  the 
V-g  method. 

A third  possible  course  is  to  apply  an  augmented-state  method, 
which  approximates  the  actual  unsteady  aerodynamic  loads  for  arbitrary 
motion  with  a transfer  function  relating  airfoil  displacements  to  loads 
having  a rational  Laplace  transform,  resulting  in  a linear  matrix  eigen- 
value problem  for  the  aeroelastic  modes.  Goland  and  Luke  (Ref.  9) 
used  this  route  to  study  wing  bending-torsion  flutter.  They  adopted 
the  R.T.  Jones  (Ref.  10)  approximation  to  the  Wagner  indicial  lift 
function  to  express  unsteady  airloads  in  rational  form,  taking  the 
Laplace  transform  in  time.  In  addition  to  their  accurate  description 
of  the  basic  bending-torsion  aeroelastic  behavior  at  all  flight  speeds, 
Goland  and  Luke  demonstrated  that  the  severity  of  flutter  cannot  be 
reliably  inferred  from  solutions  by  the  V-g  method. 
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Each  of  the  techniques  mentioned  above  attempts  to  gain  information 
about  aeroelastic  modes  whose  eigenvalues  have  nonzero  real  parts.  The 
true  effect  of  the  unsteady  wake  on  the  aerodynamic  loads  for  arbitrary 
motions  is  approximated  to  varying  degrees.  This  is  done  quite  well 
for  most  motions  in  the  case  of  augmented-state  methods,  marginally 
in  the  p-k  method,  and  not  at  all  in  the  quasisteadv  case. 

For  present  purposes  all  of  these  schemes  were  rejected  in  favor 

* 

of  the  more  exact  approach  developed  by  Edwards  (Ref.  11  ).  An  impor- 
tant contribution  of  Ref.  11  (adapted  from  Sears,  Ref.  21)  is  the 
definition  of  a generalized  Theodorsen  function  to  represent  the  exact 
circulatory  two-dimensional  incompressible  unsteady  airloads  in  the 
Laplace  domain  for  arbitrary  motions.  The  generalized  Theodorsen  func- 
tion is  expressible  in  terms  of  the  modified  Bessel  functions  of  complex 


argument 

K and 
o 

as 

Kx(s) 

(5-1) 

C(=)  E K 

(i)  + K. 

where 

sb 
s = V 


Although  previous  investigators  had  recognized  that  this  form  was 
convergent  for  the  right  half  plane,  representing  divergent  oscillatory 
motions  with  Re(s)  > 0 , the  restriction  on  the  integral  definitions 
of  the  modified  Bessel  functions  caused  some  investigators  to  believe 
that  convergent  oscillatory  motions  (s  in  the  left  half  plane)  could 

See  also  Milne  (Ref.  22). 
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not  be  so  represented.  Edwards  observed  that  K (s)  and  K, (s)  are 

o 1 

defined  and  analytic  throughout  the  s-plane  except  for  a branch  point 
at  the  origin.  When  one  places  a branch  cut  along  the  negative  real 
axis,  C(s)  can  be  shown  by  analytic  continuaLion  to  relate  circulatory 
loads  and  displacements  throughout  the  s-plane  except  along  this  cut. 

It  is,  in  effect,  an  "aerodynamic  transfer  function"  in  the  Laplace 
domain. 

With  substitution  of  s = ik  in  (5-1),  the  familiar  Theodorsen 
function  of  reduced  frequency  for  simple  harmonic  motion  is  recovered. 
Although  arbitrary  motion  is  now  being  considered  rather  than  simple 
harmonic  motion,  the  two  approaches  yield  similar  forms  when  the  initial 
conditions  arising  in  the  transforms  are  neglected.  In  fact,  the  simple 
harmonic  airloads  (4-18)  and  (2-3)  can  be  used  for  arbitrary  unsteady 
motion  simply  substituting  C(s)  for  C(k)  and  s for  icj  . 

The  modified  Bessel  functions  are  computed  from  their  ascending 
power  series  expansions,  as  mentioned  in  Ref.  11  and  described  in 
Appendix  B.  Since  the  transforms  of  aerodynamic  loads  will  be  multiple- 
valued functions  because  of  the  branch  point  of  C(s)  at  the  origin, 
the  convention 


- II  < Arg(s)  II 

is  used  for  the  cut  on  the  negative  real  axis.  The  generalized 

i0 

Theodorsen  function  is  computed  by  Edwards  in  the  form  C(fe  ) for 
representative  values  0 _<  0 _<  IT  . It  is  shown  to  approach  1 as 
r -*■  0 and  h as  r -*■  °°  for  all  0 . 
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With  the  ability  to  compute  unsteady  two-dimensional  incompressible 
airloads  for  arbitrary  motions  in  hand,  the  simple  harmonic  stability 
analysis  developed  in  Chapter  IV  can  be  generalized  for  this  case.  All 
of  the  true  aeroelastic  eigenvalues  and  eigenfunctions  can  be  obtained 
for  any  prescribed  speed  . Stability  can  be  displayed  with  root 
locus  diagrams.  The  solution  technique  is  developed  in  the  following 
section. 


B . Solution  for  the  Aeroelastic  Roots  by  Means  of  Assumed  Mo des 

Formally,  the  procedure  for  developing  the  modal  equations  needed 
for  the  true  aeroelastic  modes  throughout  the  complex  plane  begins  with 
Laplace  transformation  of  the  linearized  equations  of  motion  (3-25). 

The  transformed  perturbation  displacements  are  then  expressed  as  series 
expansions  in  the  assumed  modes  as 

n 

w (v ; s)  = £ f (y)q  (s) 

1 i-1  wi  wi 

n 

(5-2)  v.(y;s)  = Z f (v)q  (s) 

i=l  Vi  Vi 


<t>j  (v;s) 


1 f,*  <«*> 


A system  of  homogeneous,  linear,  algebraic  equations  in  the  generalized 


displacements  q 

w 

i 

method,  as  before. 


Nontrivial 


can  then  be  derived  by  Galerkin's 
solutions  are  given  by  the  zeros  of 


the 


determinant  in  s . Since  the  coefficients  in  this  determinant  which 
arise  from  the  aerodynamic  loads  will  contain  the  nonrational  function 
C(s),  this  will  not  be  a polynomial  eigenvalue  problem.  Roots 
s = a + iu>  will  thus  have  to  be  obtained  by  iteration. 
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Since  full  development  of  the  modal  equations  is  analogous  to  that 


for  the  simple  harmonic  case,  the  stability  determinant  in  s will 
simply  be  constructed  directly  from  the  modal  equations  for  simple 
harmonic  motion  defined  by  (4-12),  (4-18),  and  (2-3).  One  replaces  iu> 
by  s , ik  by  s , and  eliminates  g . The  elements  of  the  aerodynamic 
matrix  given  in  (4-18)  for  simple  harmonic  loads  become,  for  arbitrary 
motion. 


These  aerodynamic  loads  are  based  on  the  assumptions  discussed  prior 
to  Eq.  (4-13) . 

With  the  same  mass  and  stiffness  matrices  M and  K illustrated 

s s 

in  Fig.  4-3,  and  the  aerodynamic  matrix  Q whose  elements  are  defined 
by  (5-3),  the  matrix  form  of  the  modal  equations  in  s becomes 

GI. 

(5-5)  lUMj  + IQ])  + [Kg]Hq(s)}  = 0 

(5-5)  >y  be  compared  with  the  simple  harmonic  form  (4-19).  It  is 
convenient  to  define  a dimensionless  Laplace  transform  variable 


which  is  related  to  the  reduced  Laplace  transform  variable  through  the 
dimensionless  velocity  bv 

(5-7)  l = = * 

V U 

The  stability  determinant  thus  takes  the  form 
(5-8)  |p2([Msl  + [Q])  + [KJ|  = 0 

Zeros  of  this  determinant  will  yield  3n  exact  roots  for  the  aero- 
elastic  modes  in  terms  of  the  3n  assumed  modes.  These  roots  describe 
modal  frequencies  and  stability  at  the  speed  U__  used  to  calculate  the 
steady-state  deflections  which  enter  as  coefficients  in  the  stiffness 
matrix. 

A computer  program  was  developed  to  locate  numerically  the  zeros 
of  the  determinant  (5-8);  the  logic  is  outlined  in  Fig.  5-1.  This 
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format  proved  to  be  quite  convenient  for  constructing  root  locus 
diagrams  having  either  the  speed  U , the  root  angle  of  attack  a , 
or  the  drag  parameter  C as  the  changing  parameter.  Any  of  the  impor- 
tant aeroelastic  modes,  at  any  degree  of  stability,  could  be  traced 
through  the  complex  plane  as  long  as  the  initial  guess  of  s was 
sufficiently  close  to  its  particular  root  locus. 

Since  the  determinant  order  is  3n  and  n _>  2 is  desirable  to 
model  adequately  the  physical  system,  numerical  expansion  of  the 
determinant  was  not  practical.  A library  subroutine,  employed  to 
calculate  the  complex  determinant,  proved  to  be  the  source  of  numerical 
difficulties.  It  was  found  that  for  1 <_  n <_  4 and  a * 0 the 
Fig.  (5-1)  program  gave  accurate  results  when  compared  to  parallel 
V-g  method  neutral  stability  computations,  agreeing  to  at  least  seven 
digits  in  flutter  speeds.  For  n = 2 and  a t 0 , which  gives  rise  to 
steady  deflections  due  to  lift,  similar  good  agreement  was  encountered. 
But  for  n _>  3 and  a f 0 , the  program  converged  on  zeros  which 
did  not  match  the  neutrally  stable  V-g  predictions  and  were  obviously 
incorrect  from  a physical  standpoint.  Subsequent  investigation  revealed 
that  the  numerical  difficulties  originated  in  the  library  subroutine. 

When  a = 0 the  fore-and-aft  bending  degree  of  freedom  is  dynami- 
cally uncoupled,  and  the  order  of  the  determinant  which  was  actually 
computed  by  the  library  subroutine  was  reduced  to  2n  . Thus,  for 
n = 3 and  a = 0 , the  actual  computed  determinant  was  of  order  6, 
whereas  for  n = 3 and  Ot  f 0 the  order  was  9.  For  the  latter  case 
the  actual  magnitude  of  computed  determinants  was  often  0(10^)  , 
while  for  ot  = 0 the  n = 3 determinants  0(10^)  . For  n = 2 
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and  a = 0 determinants  again  were  0(10  ).  The  magnitude  of  the 

computed  determinant  thus  appears  to  be  related  to  the  numerical 
difficulties.  Since  n = 2 results  are  judged  to  model  the  problem 
adequately  and  never  encountered  numerical  problems,  correction  of  the 
above  difficulties  was  not  pursued. 

As  a result,  all  root  loci  shown  herein  for  steady  lifting  condi- 
tions involve  n = 2 . As  will  be  discussed,  however,  this  restriction 
does  not  compromise  the  modeling  of  the  physical  system  nor  prevent 
qualitative  understanding  of  its  behavior.  Furthermore,  numerical 
results  alwjiys  agree  reasonably  well  when  compared  with  V-g  computa- 
tions for  n = 3 , a t 0 . 

The  algorithm  used  to  estimate  the  zeros  of  the  determinant  in 

the  s plane,  using  an  initial  guess  s^  , is  illustrated  graphically 

in  Fig.  5-2.  The  complex  determinant  is  first  calculated  at  s^  and 

at  the  two  related  points  s + .001  and  s + i.001.  Points  A and  C 

o o 

are  then  determined,  at  which  linear  extrapolation  in  the  two  orthogonal 
directions  predicts  that  the  real  part  of  the  determinant  will  vanish. 
Similarly  points  B and  D are  predicted,  for  which  the  imaginary  part 


vanishes  by  extrapolation.  A new  guess  for  the  root  s^  is  then 
determined  as  the  intersection  point  of  dashed  lines  in  the  figure. 

The  process  is  repeated  until  satisfactory  convergence  is  realized. 
This  simple  scheme  worked  quite  well  and  never  failed  to  converge  on  a 


root,  usually  within  four  or  five  iterations.  The  convergence  criterion 
used  was  generally  ||s^|  - |s  ||  10  ^ . Typical  performance  of  the 

algorithm  is  documented  in  Table  5.1. 
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C.  Mode  Shapes  for  Aeroelastic  Modes 

The  mode  shapes  associated  with  roots  determined  for  arbitrary 
motion  by  the  determinant  iteration  method  were  conveniently  calculated 
with  the  same  linear  eigenvalue  routine  applied  in  solving  the  simple 
harmonic  stability  problem.  For  a root  computed  by  iteration,  the 
nonrational  aerodynamic  terms  containing  C(s)  in  the  matrix  0 
can  be  immediately  evaluated.  Thus  one  is  led  to  a conventional  matrix 
eigenvalue  problem,  which  obtains  the  same  root  as  one  of  its  eigen- 
values and  also  provides  its  eigenvector.  This  approach  worked  well 
and,  as  a bonus,  verified  the  accuracy  of  the  roots  computed  by  the 
determinant  iteration  scheme. 

Since  numerical  difficulties  were  encountered  with  the  determinant 
evaluation  routine  for  n = 3 and  a i 0 , this  eigenvalue  approach  of 
rechecking  its  results  for  n * 2 , a 4 0 (and  for  all  n with  ct  = 0) 
is  valuable.  It  offers  the  only  means  of  verifying  computed  roots 
lying  off  of  the  iw  axis.  Correlation  to  at  least  six  significant 
digits  was  always  observed.  Moreover,  the  accuracy  of  the  n = 2 , 
ot  # 0 determinant  iteration  solutions  has  been  checked  for  a few 
representative  cases  by  letting  n * 3 in  the  eigenvalue  routine,  with 
the  known  n = 2 root  as  a first  guess  to  evaluate  the  aerodynamic 
loads,  and  iterating  until  the  true  n * 3 root  is  obtained.  Two 
cases  for  which  this  was  done,  together  with  their  eigenvectors,  are 
shown  in  Table  5.1,  where  n = 4 roots  are  also  given.  The  n = 2 
results  obtainable  by  the  determinant  iteration  routine  are  thus  seen 
to  be  acceptably  accurate  even  for  the  higher  frequency  mode. 
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D . Inclusion  of  Unsteady  Chordwlse  Loads  Due  to  Leading  Edge  Suction 

Since  motion  of  the  wing  in  fore-and-aft  bending  is  permitted, 
unsteady  chordwise  loads  can  participate  in  the  dynamic  stability 
problem.  Unsteady  two-dimensional  incompressible  airloads  given  in 
Eqs.  (4-181  and  in  (5-3)  are  strictly  based  upon  the  assumption  that 
the  instantaneous  resultant  unsteady  lift  on  any  airfoil  section  along 
the  span  is  always  perpendicular  to  the  direction  of  the  free  stream 
velocity.  Two-dimensional  incompressible  potential  flow  theory,  however, 
does  predict  an  unsteady  leading  edge  suction  force  which  arises  from 
the  inverse  square  root  singularity  of  the  vorticitv  distribution  along 
the  airfoil  chord  at  its  leading  edge.  This  effect  will  be  included 
into  the  analysis  within  the.  framework  of  the  linearized  unsteady  per- 
turbation theory  used  to  determine  stability.  The  effect  of  the  unsteady 
propulsive  force  on  stability  can  then  easily  be  isolated  by  comparison 
of  roots  computed  for  the  (4-18)  airloads  directly  with  roots  determined 
with  the  airloads  derived  in  this  section. 

The  existence  of  a leading-edge  suction  force  due  to  the  leading- 
edge  singularity  was  determined  by  Von  Karman  and  Sears  (Ref.  12). 
Greenberg  (Ref.  13)  in  developing  the  propulsive  force  on  an  airfoil 
in  an  oscillating  stream,  states  that  a propulsive  force  acting  on  the 
airfoil  in  the  upstream  direction, 

(5-9)  F = IIoC* 

S F 


arises  from  the  unsteady  vorticitv  distribution  which  behaves  at  the 
leading  edge  as 


(5-10) 


/ b >x+ 1 


-1 
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In  this  application  the  coefficient  C will  contain  contributions 

F 

from  the  superposed  vorticity  distributions  due  to  the  steady  and  the 
unsteady  lifting  flow  fields.  Therefore  C_  will  be  a sum  of  steady 

r 

2 

and  unsteady  parts,  and  C~  will  involve  a steady  portion,  a linear 
unsteady  cross  multiplication  part,  and  a nonlinear  quadratic  unsteady 
term.  In  the  context  of  the  linearized  unsteady  perturbation  approach 
to  stability  analysis,  only  the  cross  multiplication  term  will  enter 
the  dynamic  equations.  To  include  consistently  the  nonlinear  unsteady 
propulsive  force  effect  on  stability,  the  nonlinear  structural  coupling 
terms  discarded  during  linearization  would  have  to  be  reintroduced. 

As  a consequence  of  linearization,  the  unsteady  propulsive  force 
can  be  included  only  when  both  steady  and  unsteady  vorticity  distribu- 
tions are  present.  Thus  the  case  of  zero  steady  lift  will  have  no 
contribution  due  to  this  effect  to  the  state  of  stability.  The  effect 
will  become  increasingly  pronounced  as  the  steady  lift  is  increased. 

The  vorticitv  singularity  strength  C_.  in  (5-9)  was  given  by 

1 

Garrick  (Ref.  14)  for  an  airfoil  oscillating  in  a uniform  stream  as  j 

I 

(5-11)  C = /2b  { [h  + Vex  + ba('i-a)  ]C(k)  - ’jba} 

r 

with  h positive  downward.  Converting  to  present  notation,  introducing 
the  generalized  Theodorsen  function,  and  introducing  superposition  of 
steady  and  unsteady  deflections  gives 

a -*■  a + <(>o  + <J>^ 

h ■*  - w - w, 
o 1 
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(5-12)  C = /2b  { (a+(J>  )V  + [-w  + <)>  V + b^-a)^.  ]C(s)  - ’&<{)..} 

r Oil  i 1 


The  propulsive  force  Fg  , using  (5-9),  is 

(5-13)  F = 2Jlpb{  (a+4>  )2  V2  + 2V(a-HJi  )[(-/  + <J>,V 
s o oil 

+ b(b-a)<51)C(s)  - '.bc^]  + [ (-w^  - <^V  + b(?5-a)({)1)C(s) 

- ,ib<j>1  ]2 

The  last  squared  term  is  the  nonlinear  time-dependent  contribution  and 
is  neglected  hereafter.  The  first  squared  term  is  the  propulsive  force 
on  a flat-plate  airfoil  at  incidence  in  steady  flow.  The  resultant 
steady  aerodynamic  force  should  act  at  right  angles  to  the  free  stream 
velocity  in  potential  flow,  and  this  steady  propulsive  force  can  be 
interpreted  physically  as  the  component  which  tilts  the  resultant  lift 
vector,  obtained  by  summing  the  pressure  distribution  at  right  angles  to 
the  chord,  forward  to  become  normal  to  the  airstream. 

The  propulsive  force  can  thus  be  seen  to  correct  for  the  chordwise 
component  of  the  lift  which  is  computed  normal  to  the  airfoil  chord. 

The  assumption  incorporated  into  (4-18)  that  the  unsteady  circulatory 
lift  acts  at  right  angles  to  the  airstream  must  be  discarded  and  the 
force  assumed  rather  to  be  normal  to  the  chord  of  the  airfoil  in  its 
steady-state  orientation.  The  assumed  direction  of  forces  is  shown  in 
Fig.  5-3. 

In  the  structure  axis  system,  the  loads  are 
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(5-14) 


F cosd>  + L slnd) 
s o a o 


F =1,  eosd>  + F si  nd) 
z . a o s i> 


flu*  following  assumptions  arc  made 


COSif)  ' I 


s imj)  ~ A 
o o 


F si  nit)  <<  L 

s o 


so  that 

F “ L 
71  a 

F = - F + I , i]i 

x|  8 8 0 

The  lift  for  arbitrary  motion  in  lime  is  given  bv  Kef.  (ri)  as 
<5-16)  b;)  = 2HpbV[-  w(  C(s)  + (J)(  V C(s)  + bCj-a)^!  C(s)  ] 

+ llpb2f-w1  + - ba<f)  1 

Combining  (r)-l(>)  and  the  linear  unsteady  part  or  (5-13)  into  the  second 
of  (5-15)  yields 


- -4IIpVb (<x  + ’.4,)  I-  wf  C(s)  + 1 V C(s) 


+ ((Vn)C(s) 


’s)b<b,  I + llpirvi})  <j> 

I O 1 


+ npb'iM-w.  + Vif*  - baif>.  I 


(5-17) 


The  loads  can  now  be  arranged  into  the  unit  generalized  force 
format  of  (5-3)  hv  Laplace  transformation  in  time,  substitution  of  the 
assumed  modes  as  in  (5-2),  and  formation  of  the  generalized  forces  from 
relations  such  as  (4-11),  After  this  work,  which  is  straightforward, 
*-he  elements  of  the  aerodynamic  matrix  which  incorporate  the  linearized 
unsteady  propulsive  force  are 


(5-18) 


V,i  * 2<1-Va  V, 


n 

+ (2-L  ) T.  H , q" 
w . , vx 
V=m+ 1 J 


Q.,  = 2[L.  - - •=■  + A(l-L  ) ]I,.  a 

j+n,i+2n  <p  s w ji 


(1<  i<  n) 
(1<  j<n) 


+ [L.  - 1 - + (2-L  ) A 1 Z Y . q° 

4>  s w ,,  Vi  V 

b V=m+1  J v 


where  the  new  modal  integrals  appear. 


(5-19) 


l 

Y = / f f.  ft 
ivj  o W 


dv 


The  remainder  of  the  terms  of  the  aerodynamic  matrix  remain  the  same 
as  in  (5-3). 


Prior  to  actual  calculations  a further  approximation  is  made. 


The  terms  in  (5-18)  which  depend  on  sums  of  q 


are  neglected. 


eliminating  the  need  to  compute  the  . This  is  equivalent  to 

assuming  that  the  lift  L in  Fig.  5-3  is  aligned  with  the  z-axls  and 

(1 

that  <{>^  is  neglected  relative  to  a in  the  linear  term  of  (5-13). 
Strictly,  this  simplification  will  alter  the  results  somewhat,  but  it 


is  not  expected  to  change  the  overall  effect  of  the  propulsive  force 


on  stability  and  does  simplify  computations.  The  first  order  trend  of 
the  effect  of  100%  leading  edge  suction  on  stability  is  the  main  point 
of  interest  and  should  not  be  affected. 

The  changes  made  to  the  aerodynamic  matrix  are  therefore  substitu- 
tion of  the  terms 


(5-20) 


Q 


j+n.i 


< 2(1-L  )a  (i=j) 

w 

0 <i*J> 


Q.,  = 2 [L,  - h - - + A(l-L  )]I  a 

j+n* i+2n  1 0 s w ji 


(1<  j <n) 
(1<  i<n) 


for  their  counterparts  in  (5-3) . 

The  program  described  in  Section  B includes  the  option  of  using 
either  of  these  unsteady  aerodynamic  force  systems,  and  a comparison  of 

their  relative  effect  on  stability  is  made  in  the  next  chapter.  Except 

2 

for  the  - — term  in  the  second  of  (5-20),  incidentally,  the  newer 

s 

system  simply  involves  doubling  the  magnitude  of  the  terms  in  (5-3) 
that  are  replaced. 
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M = 40,  P = 0.005,  ia  = 0.25,  A = 0.1,  S = 0.1,  T = 25 , a = 0 , n = 3, 
C = 0 , U " 6 . 5 


ITERATION 

s 

DETERMINANT  AT  s 

0 

-0.072968  + iO. 0582071 

0.326xl09  + i0.1266xl010 

1 

-0.078290  + 10. 067052 

-0.259xl09  + 10.504xl08 

2 

-0.079541  + 10. 065839 

0.106xl08  + i0.504xl07 

3 

-0.079527  + iO. 065910 

0.106xl06  - i0.178xl06 

4 

-0.079526  + iO. 065910 

-0.281xl04  + i0.499xl04 

M = 40,  P = 0.02,  in  *=  0.25,  A = 0.1,  S -0.1,  T = 25,  a = 0.02,  n = 2, 
C = 0,  U = 7 


ITERATION 

s 

DETERMINANT  AT  s 

0 

-0.0055179  + iO. 119612 

0.261xl0U  - i0.253xl010 

1 

-0.0014167  + iO. 1074866 

0.125xl010  + i0.428xl010 

2 

0.0005501  + iO. 1085285 

-0.459xl07  + I0.135xl09 

3 

0.0005027  + iO. 1084861 

-0.129xl06  + i0.178x!06 

4 

0.0005026  + iO. 1084861 

TABLE  5.1  Performance  of  the  Determinant  Iteration  Algorithm  for  Two 
Cases,  One  Nonlifting  With  n = 3 and  One  at  Steady  Lift 
With  a = 0.02,  n * 2 
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STEADY- ST AT E DEFLECTION S : 


C /b 
q:./b 
<C  /b 

q°3/b 

w 

4 

< /b 
<1  ° 1 /b 
V2 

q°"/b 

V3 

qol/b 

o4 

qti 

’*7 


n = 2 
1.97768 
0.022566 


0.003170 

-0.00363 


0.0068192 

-0.0001112 


n = 3 
2.01496 
0.023221 
0.001453 

0.003534 

-0.000403 

-0.000044 

0.0072034 

0.0000129 

-0.0000322 


n = 4 
2.019088 
0.023274 
0.001464 
0.000271 
0.003594 
-0.000407 
-0.000046 
-0.000007 
0.0072456 
0.0000303 
-0.0000130 
-0.0000103 


M - 40 

P * 0.005 

1 =0.25 

a 

A = 0.1 
S =0.1 
T = 60 

a = 0.01  RAD. 
C = 0 
U = 7 


ROOT  FOUND  BY  DETERMINANT 

ITERATION  WITH  n = 2:  p = 0.03579  + 10.55875 

RESULTS  OF  LINEAR  EIGENVALUE  ANALYSIS  FOR  n = 2,  3,  4: 


n = 2:  n = 3: 

p = .035795  + 1.55875  p = .036828  + 1.52903 

]-  T 

AMPLITUDE1  PHASE  AMPLITUDE1  PHASE 


v/b 

qw  /b 

- 2 

q /b 

-3 
q /b 
w 

q 4 /b 

-1 
q /b 

>2 

q /b 

-3 

V” 

>3 


4.2245 

1.2064 


0.083839 

0.05128 


1.0 

0.05780 


PHASE 

AMPLITUDE  1 
| 

PHASE 

215.40° 

4.6573  1 

217.75 

1 

o 

o 

1.1808  ' 

i 

-37.53 

0.03432 

-18.25 

- 7.71° 

0.87304  | 

- 8.79 

181.23° 

0.05146  | 

181.49 

0.009128  1 
| 

169.96 

0.0° 

1.0  1 

I 

0.0° 

-33.85° 

0.07738 

-30. 30 

0.02511  | 

1 

-19.79 

n = 4 : 

p = .036486  + 1,52309 
AMPLITUDE!  PHASE 


4.7497 

1.1815 

0.03544 

0.00530 

0.88887 

0.050487 

0.009314 

0.002135 

1.0 

0.08213 

0.02969 

0.00912 


218.17° 
-37.07° 
-17.97° 
-12.46° 
- 9.12° 
181.93° 
169.12° 
165.92° 
0.0° 
-29.93° 
-18.96° 
-22.80° 


TABLE  5.2  Modal  Convergence  for  Convergent  and  Divergent  Oscillatory 
Aeroelastic  Modes  (see  Fig.  6-24(a)) 
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ROOT  FOUND  BY  DETERMINANT 

ITERATION  WITH  n = 2:  p = -0.216505  + 11.632293 


RESULTS  OF  LINEAR  EIGENVALUE  ANALYSIS  FOR  n = 2,  3,  4 


n = 2: 

n = 3* 

n = 4 : 

p = -.2165053  + 

11. 6322931 

p = -.21485  + 

11.59221 

p = -.21564  + 

11.57794 

| 

AMPLITUDE  1 

1_ 

PHASE 

AMPLITUDE ' 

PHASE 

~T 

AMPLITUDE | 

PHASE 

q /b 

0.02478  1 

175.84° 

0.02627 

175.18° 

0.025894 1 

174.95° 

1 

1 

1 

q / b 

w 

2 

0.05039  | 

181.31° 

0.07554  | 

1 

178.03° 

0.07883  , 

177.63° 

q /b 
w 

3 

0.02535  1 

1 

181.07° 

0.02054  | 

1 

178.99° 

q /b 
W4 

1 

0.00047  1 

159.19° 

qv  /b 
_ l 

qv  /b 

0.14780  1 

| 

160.48° 

0.14548  | 

160.42° 

0. 14364 

160.30° 

0.04973  , 

181.88° 

0.05680  | 

181.47° 

0.05937  | 

181  . 34° 

V2 

i 

q /b 

| 

0.000618  1 

122.60° 

0.000565  1 

242.35° 

V3 

1 

1 

t 

q /b 

1 

1 

0.000081 

- 6.87° 

V4 

| 

1 

1 

\ 

1.0 

0.0° 

1.0  I 

1 

0.0° 

1.0 

0.0° 

0.20484  | 

175.75° 

0.18586  1 

1 

175.28° 

0.18223  | 

175.24° 

I 

0.05215  j 

176.65° 

0.05252  1 

1 

176.50° 

\ 

1 

L 

L 

0.01955 

1 

177.66° 

TABLE  5.2  CONCLUDED 
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Chapter  VI 

AEROELASTIC  MODES  USING  AIRLOADS  FROM  INCOMPRESSIBLE  STRIP  THEORY 


A.  The  Effect  of  Steady  Drag  on  Flutter  of  a Nonlifting  Wing 

Before  considering  steady-state  deflections  due  to  lift,  a thorough 

understanding  of  the  stability  behavior  of  the  cantilever  wing  at  zero 

steady  lift  is  needed.  With  w and  d>  at  zero,  the  fore-and-aft 

o o 

bending  degree  of  freedom  is  dynamically  uncoupled  from  vertical  bending 
and  torsion  motions,  and  the  system  analyzed  in  Ref.  1 results.  In  the 
zero  lift  case,  then,  solutions  for  stability  involve  2n  aeroelastic 
modes  consisting  of  coupled  motions  in  and  <j>^  ; the  remaining  n 

modes  represent  uncoupled  free  vibration  in  each  of  the  assumed  modes 
in  v ^ . As  demonstrated  in  previous  chapters,  the  flutter  conditions 
obtained  by  this  assumed  mode  analysis  compare  favorably  with  Ref.  1 
results  over  all  practical  combinations  of  the  parameters  M , P , i , 

A , S , and  C . Owing  to  this  good  agreement,  the  results  and  conclu- 
sions of  Ref.  1 apply  here  as  well,  yet  the  assumed  mode  solution  method 
still  is  useful  in  providing  additional  insight  into  the  flutter  behavior 
of  the  nonlifting  wing. 

The  parameters  M , i , A , and  S offer  no  suprising  effect, 
and  most  importantly  here,  the  nature  of  their  influence  is  not  altered 
by  the  inclusion  of  drag.  The  Ref.  1 results  indicate  that  an  increase 
in  the  elastic  axis  - A.C.  offset  given  by  A is  destablizing,  an 
increase  in  the  sectional  C.G  - elastic  axis  offset  given  by  S is 
destabilizing,  an  increase  in  the  radius  of  gyration  parameter  i is 
stabilizing,  and  that  the  flutter  speed  is  approximately  proportional 
to  the  square  root  of  the  mass  ratio  parameter  M . Since  Ref.  1 
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establishes  that  the  influence  of  drag  is  not  sensitive  to  these  four 
parameters,  more  detailed  study  in  this  area  is  not  considered  here. 

The  most  curious  finding  of  Ref.  1 involves  the  effect  of  C in 
conjunction  with  the  aspect-ratio  parameter  P , which  is  the  product 
of  the  vertical  bending-to-torsion  stiffness  ratio  with  the  inverse 
square  of  the  geometric  aspect  ratio.  The  effect  of  steady  drag  on 
flutter  speed  is  stabilizing  for  smaller  aspect  ratios  (larger  P)  but 
destabilizing  for  larger  aspect  ratio  wings.  The  reversal  of  the  effect 
of  drag  on  flutter  occurs  near  P - 0.01  , and  the  behavior  in  this 
neighborhood,  including  flutter  mode  shapes,  appears  to  be  quite 
interesting.  No  conclusions  regarding  physical  causes  of  this  phenomenon 
were  made  in  Ref.  1,  however. 

Figure  6-1  is  a reproduction  from  Ref.  1,  showing  the  effect  of 
drag  on  flutter  speed  as  a function  of  P for  intermediate  values  of 
M , i^  , A , and  S . Numbers  in  parentheses  on  the  abscissa  give  the 
true  aspect  ratio  for  a typical  value  of  the  ratio  EI^/GI^  = 1.6  . 

Clearly  wings  of  practical  interest  include  the  region  within  which  the 
effect  of  drag  on  flutter  appears  to  be  most  interesting. 

To  help  gain  a better  physical  understanding  of  the  behavior  near 
P = 0.01  , flutter  solutions  for  this  same  example  have  been  found  via 
the  simple  harmonic  method  of  Chapter  IV,  over  the  range  0.002  <_  P 0.02. 
Results  appear  in  Fig.  6-2  that  show  flutter  speeds  and  flutter  mode 
shape  amplitudes  and  phase  relations  as  functions  of  P . Since  a 
finite  chordwise  to  vertical  bending  stiffness  ratio  T must  be 
specified,  and  the  effect  of  drag  on  flutter  depends  on  T as  in 
Fig.  4-15  the  value  T * 50  was  used  to  allow  adequately  for  the 


T -*■  00  behavior  inherent  in  the  Ref.  1 formulation.  Three  assumed 
modes  in  each  degree  of  freedom  are  used. 

The  flutter  speeds  in  Fig.  6-2a  closely  match  the  Fig.  6-1  results. 

The  flutter  mode  shapes  include  participation  by  the  three  generalized 

displacements  q,  , q , and  q , with  the  remaining  assumed  modes 

*1  W1  W2 

contributing  negligibly  to  the  motions.  In  Fig.  6-2b  the  amplitude  and 
phase  of  the  two  assumed  bending  modes  at  flutter  are  shown  for  unit 
magnitude  and  zero  phase  angle  of  the  first  assumed  torsion  mode. 

It  can  be  seen  that  for  any  C the  sharp  drop  in  flutter  speed 
that  comes  with  decreasing  P is  accompanied  by  a sudden  change  in  the 
flutter  mode  shape.  The  amplitudes  of  the  two  bending  modes  merge,  and 
the  second  bending  mode  undergoes  a large  phase  shift.  Further  decrease 
in  P gives  a gradual  separation  of  the  assumed  bending  mode  amplitudes 
with  the  fundamental  mode  again  becoming  dominant. 

To  help  visualize  the  physical  appearance  of  these  flutter  mode 
shapes,  phasor  diagrams  of  the  spanwise  distribution  of  bending  displace- 
ments, for  C = 0.02  , are  given  in  Fig.  6-3  for  five  values  of  P . 
Arrows  depict  the  modal  generalized  displacements  from  Fig.  6-2,  and 
the  curves  give  the  relative  displacements  along  the  span  and  their 
phase  referenced  to  q^  . For  both  P = 0.02  and  P = 0.002  all 
stations  are  nearly  in  phase  and  the  mode  shape  is  dominated  by  the 
first  assumed  bending  mode.  But  for  the  intermediate  values  of  P , 
where  the  transition  in  phase  of  the  second  bending  mode  takes  place, 
the  displacements  at  different  locations  along  the  span  can  be  over  90° 


out  of  phase. 
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The  behavior  of  flutter  modes  in  this  range  of  P offers  the 


greatest  discrepancy  found  between  the  results  reported  in  Ref.  1 and 
the  assumed  mode  solution,  suggested  by  Fig.  2-1.  For  P < 0.04  in 
Ref.  1,  collocation  at  only  five  spanwise  stations  was  used,  and  the 
mode  shapes  in  bending  and  torsion  were  permitted  to  have  spanwise  phase 
differences.  It  was  found,  however,  that  these  phase  differences  never 
exceeded  a few  degrees,  in  contrast  to  the  results  presented  here. 
Possibly  the  use  of  only  five  spanwise  collocation  points  did  not  allow 
enough  freedom  to  represent  the  flutter  mode  shape  transition  found 
using  assumed  modes.  In  any  case  good  agreement  between  flutter  speeds 
and  frequencies  is  still  observed  for  the  two  methods. 

Figure  6-2  seems  to  indicate  that  the  second  assumed  bending  mode 
plays  a significant  role  in  the  reversal  of  the  effect  of  drag  on  flutter 
near  p = .01  , which  coincides  with  the  natural  frequency  of  this  mode 
crossing  the  flutter  frequency.  Interaction  of  actual  aeroelastic  modes 
is  masked  by  the  limitations  of  the  solution  method,  however,  which  only 
gives  neutrally  stable  solutions.  In  order  to  better  understand  these 
results,  the  Laplace  transform  approach  of  Chapter  V is  used  to  allow 
tracing  all  of  the  aeroelastic  modes  in  the  complex  plane  for  speeds 
from  zero  into  the  supercritical  range. 

Figures  6-4 (a)  - 6-4 (h)  give  root  locus  diagrams  for  increasing 
speed  at  eight  representative  values  of  aspect-ratio  parameter  P . 

F.ach  locus  originates  for  U = 0 at  one  of  the  normal  modes  of  free 
vibration  of  the  structure,  which  are  easily  calculated  in  terms  of  the 
uncoupled  assumed  bending  and  torsion  modes.  Zero  drag  branches  are 
shown  in  all  of  the  figures,  with  loci  corresponding  to  C ^ 0 added 
where  their  behavior  differs  significantly  from  that  for  zero  drag. 
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The  first  two  figures  represent  stubby  low-aspect-ratio  wings,  for 
which  the  strip-theory  aerodynamic  assumption  is  certainly  inadequate. 

Due  to  the  manner  in  which  the  Laplace  variable  p is  nondimensionalized , 
the  predominantly  first  torsion  normal  mode  of  free  vibration  remains 
essentially  fixed  on  the  iu>  axis  near  1.6  on  all  of  these  diagrams. 

The  predominantly  bending  normal  modes  move  down  the  iu)  axis  as  P 
decreases,  since  their  natural  frequencies  decrease  relative  to  the 
torsion  frequencies.  In  Fig.  6-4(a)  the  bending  branch  of  the  aero- 
elastic  modes  leads  to  flutter,  whereas  in  Fig.  6-4(b)  the  torsion 
branch  eventually  becomes  unstable.  The  normal  mode  having  the  third 
lowest  natural  frequency,  predominantly  the  second  bending  assumed  mode, 
occurs  well  up  the  iu>  axis  and  off  these  two  diagrams  and  has  negligible 
influence  on  flutter.  These  low-aspect-ratio  cases  show  entirely  two- 
degree-of-freedom  behavior  and  closely  resemble  the  root  locus  given 
by  Edwards  (Ref.  11)  for  a typical  section  in  plunge-and-pitch  motion 
in  incompressible  flow. 

Figures  6-4 (c)  and  6-4 (d)  represent  values  of  P just  above  the 
condition  where  the  effect  of  drag  on  flutter  reverses.  Although  the 
flutter  phenomenon  is  still  similar  to  that  for  larger  P , the  third 
normal  mode  frequency  has  now  decreased  sufficiently  to  appear  on  the 
diagram,  and  it  produces  a branch  which  does  not  lead  to  flutter  for  all 
values  of  drag. 

According  to  the  flutter  curve  of  Fig.  6-2,  for  P = 0.01  the 
cases  C = 0.02  and  C = 0.04  result  in  decrease  of  flutter  speed 
from  the  zero-drag  condition  whereas  for  C = 0.01  it  is  still  increased. 
Figure  6-4 (e)  gives  the  P = 0.01  root  locus,  which  reveals  that  the 
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aeroelastic  mode  emanating  from  the  predominantly  second  bending  normal 
mode  now  plays  an  important  role.  The  next  illustration.  Fig.  6—4 ( f ) 
with  P = 0.007,  gives  this  branch  as  becoming  unstable  for  all  C . 

In  the  final  two  of  these  illustrations,  with  P = 0.005  and  0.002  , 
the  flutter  phenomenon  appears  to  be  returning  to  the  type  of  behavior 
seen  for  small-aspect-ratio  wings,  with  the  second  bending  contribution 
assuming  a lesser  influence.  In  Fig.  6-4(h)  a fourth  normal  mode,  the 
third  bending  mode,  has  made  its  appearance  but  does  not  noticeably 
influence  flutter. 

The  nature  of  the  aeroelastic  modes  on  the  various  branches  of  the 
root  locus  diagrams  can  be  clarified  by  looking  at  their  mode  shapes. 

In  Fig.  6-5,  phasor  diagrams  are  used  to  show  C = 0 mode  shapes  for 
each  branch  of  the  .005  locus  (Fig.  6-4(g)).  At  selected  speeds  both 
subcritical  and  supercritical,  the  generalized  displacements  are  shown 

with  the  phase  angle  of  q taken  as  zero. 

<?1 

The  dependence  of  drag's  effect  on  flutter  upon  P is  also  given 
for  a second  configuration  having  a smaller  mass  ratio  M = 9.4  , more 
representative  of  light  aircraft  and  sailplanes.  Shown  in  Fig.  6-6, 
the  basic  behavior  resembles  the  first  configuration,  with  certain 
differences.  For  example,  the  reduction  in  flutter  speed  by  drag  for 
small  P is  moderate  relative  to  that  for  the  larger  mass  ratio;  this 
is  also  seen  in  the  tabulated  Ref.  1 results.  Also  the  transition  as 
P reduces  through  the  0.01  region  is  much  less  severe.  Again  the 
reversal  coincides  with  the  crossing  of  the  second  assumed  mode  frequency 
below  the  flutter  frequency. 

In  conclusion,  assumed  mode  results  attribute  the  reversal  of  the 
effect  of  drag  on  flutter  speed  to  the  interaction  of  the  structure's 

J 
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second  natural  bending  mode  with  t ho  basic  cantilever  wing  flutter 
mechanism  involving  the  first  bending  and  first  torsion  modes.  The 
exact  reason  that  drag  increases  the  flutter  speed  for  1’  above  this 
apparent  "resonance"  condition  and  decreases  It  for  smaller  1’  Is  not 
evident  in  these  results,  and  like  many  aetoelastic  phenomena  may  not 
have  a simple  physical  explanation.  It  does  appear,  however,  that  the 
drag  force  enhances  the  coupling  ol  the  second  bending  mode  Into  the 
flutter  behavior  and  thereby  magnifies  Its  already  present  effect  on 
flutter  speed. 

Since  the  Kef.  I tabulated  results  give  a fairly  complete  picture 
ol  flutter  ol  the  nonllftlng  case  In  the  presence  of  steady  drag  and 
are  not  disputed  by  current  results,  further  work  herein  is  directed 
towards  the  more  general  case  involving  steady  deformations  due  to  lift. 

II.  Effect  ol  Steady  Deformat  ions  Due  to  l.il  t 

The  good  agreement  in  prediction  of  dynamic  stability  between  the 
current  analysis  and  the  Ref.  I collocation  method  for  nonlifting  wings 
with  a steady  drag  force  included  furnishes  confidence  that  the  modal 
scheme  will  be  successful  for  steady  lift ing  conditions.  The  effect  ol 
steady  deformations  due  to  lift  is  incorporated  into  the  dynamic 
stability  analysis  by  the  same  means  as  the  steady  drag  effect — namely 
through  coefficients  ol  the  stiffness  matrix  determined  in  a separate 
nonlinear  solution  for  the  steady-state  defied  ions.  Thus,  the  agreement 
indicates  that  the  scheme  allowing  small  lime  dependent  perturbations 
about  a static  deflection  is  working  properly.  Wien  steady  lifting 
deflections  are  introduced  through  a nonzero  root  angle  ol  attack  a , 


the  fore-and-aft  bending  is  no  longer  dynamically  uncoupled. 

Coupling  in  both  elastic  terms  (the  stiffness  matrix  coefficients 


containing  and  and  aerodynamic  coupling  terms  (arising  from 

unsteady  aerodynamic  force  components  in  the  x-direction)  now  appears 
in  the  fore-and-aft  bending  dynamic  equations.  The  3n  aeroelastic 
modes  will  consist  of  coupled  motions  in  , (j^  , and  . 

Three  basic  wing  configurations  are  selected  to  illustrate  the 
effects  introduced  by  steady-state  lifting  deformations.  Parameters 
M , i , A , and  S are  taken  the  same  as  in  Fig.  6-1,  and  aspect  ratio 
parameter  P is  assigned  three  different  values  in  order  to  consider 
wings  of  large,  moderate,  and  low  aspect  ratio.  For  large  aspect  ratio, 
P = 0.005  is  chosen  to  provide  a case  for  which  steady  drag  decreases 
the  flutter  speed  (Fig.  b-1).  A moderate  aspect  ratio  example  with 
P = 0.02  having  an  increase  in  flutter  speed  due  to  steady  drag,  and  a 
low-aspect-ratio  case,  P = 0.1  , are  also  included.  For  the  typical 
bend ing-to-torsion  stiffness  ratio 


^T1'6 


these  examples  correspond  to  actual  retangular  planforms  having  aspect 
ratios  of  17.89  , 8.94  , and  4 respectively.  The  bending  stiffness 
ratio  T now  becomes  important,  and  is  given  the  nominal  value  25. 

The  essential  features  of  the  flutter  behavior  encountered  when  steady 
deformations  enter  can  be  illustrated  by  using  these  basic  configurations 
as  examples. 

Tlie  flutter  stability  boundary  for  the  moderate-aspect-ratio  example 


is  shown  in  Fig.  6-7  for  C = 0 and  C = 0.01;  the  steady  bending 
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displacement  w^  of  the  wingtip  is  the  measure  of  steady  lift.  An 

alternative  would  be  to  show  flutter  speeds  as  a function  of  a , but 

this  Is  a poor  means  for  comparing  curves  having  different  steady  drag 

and  gives  no  information  about  the  elastic  steady  deformations.  A 

better  way  to  indicate  the  steady  flight  condition  would  be  the  total 

lift  force  on  the  deformed  wing  nondimensionalized,  for  example  as 
L<1? 

• As  can  be  seen  in  Fig.  6-8  this  dimensionless  total  lift  para- 
x 

meter,  which  depends  on  the  steady  twist  <^(y)  . varies  for  constant 
wingtip  deflection  as  C changes  for  points  along  the  stability 
boundary  of  Fig.  6-7.  This  is  because  the  drag  force  alters  the  rela- 
tive i|>  and  w distributions  for  the  same  total  steady  lift.  But 
o o 

since  this  effect  is  small,  and  w gives  the  best  indication  of  the 

o 

magnitude  of  the  steady  equilibrium  deflections,  this  deflection  is  used 
to  Indicate  the  steady  lift  condition.  In  Fig.  6-7,  the  semispan  of  the 
wing  is  about  9 semichords,  and  steady  deformations  well  exceeding  the 
limits  of  the  moderate  displacement  beam  theory  are  therefore  shown. 

This  demonstrates  that  flutter  solutions  can  be  found  for  arbitrarily 
large  steady  displacements  and  that  it  is  a matter  of  practical  engineer- 
ing Judgment  to  recognize  when  the  assumptions  made  in  the  derivation 
of  the  equations  have  been  violated. 

For  C = 0 Fig.  6-7  indicates  that  the  flutter  speed  reduces 
continuously  with  increasing  steady  lift  until  a maximum  reduction  of 
about  13%  is  achieved  at  an  excessively  large  deflection  of  4 semichords. 
The  frequency  at  flutter  reduces  monotonically  with  w^  ; this  effect 
is  generally  observed  for  all  wing  configurations. 
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With  steady  drag  included,  the  situation  for  the  deformed  wing  is 


not  as  simple  as  it  is  for  the  a - 0 case.  As  exhibited  by  the  non- 
linear steady  deformations  shown  in  Chapter  IV,  drag  can  greatly  alter 
the  deformation  state  associated  with  a given  speed  and  root  angle  of 
attack.  Even  more  importantly,  drag  significantly  reduces  divergence 
speeds,  possibly  to  less  than  the  flutter  speed.  This  is  revealed  in 
the  stability  analysis  when  the  nonlinear  steady  solution  blows  up 
before  dynamic  instabilities  appear. 

In  Fig.  6-7  divergence  speeds  found  by  the  linear  VT  determinant 
(2-29)  are  also  indicated  for  several  values  of  C . Results  of  a 
dynamic  stability  analysis  with  C ■ 0.01  also  appear.  The  flutter 
speed  for  a = 0 and  C “ 0.01  is  the  single  point  on  the  ordinate 
and  is  greater  than  the  divergence  speed  for  the  same  drag.  During  the 
search  for  a neutrally  stable  oscillating  condition  for  the  very  small 
angle  of  attack  a « 0.001  rad.  and  C ■ 0.01  , the  nonlinear  steady 
displacement  solution  was  sensitive  to  U in  the  neighborhood  of  flutter 
owing  to  the  proximity  of  divergence.  For  specified  larger  angles  of 
attack,  the  C » 0.01  flutter  boundary  is  found  without  difficulty. 

The  slight  increase  in  flutter  speed  due  to  drag  observed  for  a - 0 
appears  to  be  preserved  in  the  presence  of  steady  deformations.  For 
decreasing  ot  , however,  a point  is  reached  as  a 0 at  which  the 
steady  displacements  are  still  nonzero  at  flutter.  For  smaller  steady 
deflections  divergence  becomes  the  instability  encountered  for  increasing 
speed.  For  C = 0.02,  divergence  was  observed  for  all  lifting  conditions 
with  no  flutter  points  found. 
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As  steady  deformation  increases  from  zero,  the  flutter  mode  shape 
undergoes  a smooth  transition,  originating  with  the  same  zero  lift  mode 
shapes  shown  in  Fig.  6-2  for  P = .02  . For  zero  drag,  the  amplitudes 
of  modal  generalized  displacements  which  contribute  significantly  to 
the  flutter  mode  are  shown  in  Fig.  6-9,  normalized  to  the  first  torsion 
mode  amplitude.  Their  phase  angles  relative  to  zero  phase  for  q 


appear  in  Fig.  6-10.  Participation  by  the  first  chordwise  bending  mode 
increases  steadily  with  increasing  w^  , reflecting  the  increased  strength 
of  the  elastic  bending-torsion  coupling.  Vertical  bending  motions  are 
increasingly  dominated  by  the  first  assumed  bending  mode,  and  the 
contribution  of  the  second  torsion  assumed  mode  increases  significantly. 

The  flutter  mode  shapes  at  one  steady  lifting  condition  are  presented 
in  Fig.  6-11  in  a form  giving  a clearer  physical  description  of  the 
motion.  Above  the  phasor  diagram  (which  contains  the  information  given 
in  Fig.  6-9  and  6-10)  is  a sketch  of  the  cyclic  path  traced  in  the 
x-z  plane  by  wing  sections  at  the  wingtip  and  at  midsemispan.  Points 
where  the  first  torsion  assumed  mode  is  at  phase  angles  of  0°  , 90°  , 

180°  , and  270°  are  located.  These  diagrams  emphasize  the  three-degree- 
of-freedom  nature  which  flutter  can  have  when  steady  deflections  are 
present . 

The  low-aspect-ratio  example  (P  = 0.1  and  /$?.  ~ 4)  shows  only  minor 
effects  upon  its  flutter  characteristics  due  to  steady  deformations,  as 
might  be  anticipated.  As  given  in  Fig.  6-12,  even  for  the  extreme 
condition  a = . 12  rad.  yielding  a 1.7  semichord  tip  deflection  at 
flutter,  there  is  only  a 0.62%  reduction  in  flutter  speed  due  to  lift. 
The  flutter  mode  shapes  (Fig.  6-12)  undergo  little  change,  with  a slight 
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contribution  by  q the  only  new  feature.  Owing  to  these  unremarkable 
V1 

results,  further  work  centers  on  the  moderate-and-high-aspect-ratlo 
examples,  where  the  effects  of  steady  deformations  and  chordwise  forces 
are  significant  and  interesting. 

For  the  high-aspect-ratio  example  with  P = 0.005  , the  flutter 
behavior  is  quite  different  when  steady  deformations  enter.  Figure 
6-13  shows  the  dependence  of  flutter  speed  upon  the  steady  wingtip 
bending  deflection  for  two  drag  cases,  C « 0 and  0.01  . The  minimum 
flutter  speed  in  this  case  is  over  20%  below  its  undeflected  counterpart, 
but  quite  interestingly  as  w ^ 0 the  stability  boundaries  do  not 

converge  continuously  to  their  respective  zero  lift  flutter  speeds  but 
approach  lower  points  on  the  ordinate. 

The  flutter  mode  shapes.  Figs.  6-14,  6-15,  explain  this  new 
behavior.  As  the  steady  deflection  becomes  small,  the  flutter  mode 
becomes  dominated  by  the  first  chordwise  bending  mode,  and  as  a -*■  0 
titis  type  of  instability  approaches  simple  free  uncoupled  vibration  in 
this  degree  of  freedom.  For  moderately  large  steady  deflection  the 
flutter  mode  shape  closely  resembles  that  for  the  medium-aspect-ratio 
wing.  Figs.  6-9  and  6-10. 

Flutter  mode  shapes  for  two  steadv-lift  conditions  are  diagrammed 
in  Fig.  6-16  using  the  same  technique  as  in  Fig.  6-11.  Relative  to  the 
P = 0.02  example,  this  wing  shows  a greater  amount  of  participation  in 
fore-and-aft  bending,  and  a greater  contribution  from  the  second  vertical 
bending  mode.  This  mode's  contribution  causes  the  noticeable  difference 
in  the  eccentricities  of  the  elliptical  paths  traced  in  the  x-z  plane 


bv  different  wing  stations.  The  phase  relationship  among 
and  q appears  similar  in  this  and  other  examples. 

v . 
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The  most  noticeable  effect  of  steady  drag  again  is  its  reduction 
of  the  divergence  speed.  Figure  6-17  gives  divergence  speeds  found  by 
the  linear  determinant  (2-29)  together  with  zero  lift  flutter  speeds 
for  increasing  steady  drag  for  the  moderate-  and  high-aspect-ratio 
cases.  Divergence  is  clearly  more  important  relative  to  flutter  for 
the  larger  aspect  ratio. 

The  effect  of  steady  drag  on  the  dynamic  behavior  (Fig.  6-13) 
appears  to  cause  only  slight  adjustments  to  the  zero  drag  flutter  results 
for  any  steady  deflection.  The  flutter  mode  shape  amplitudes  and  phase 
angles  show  but  small  changes  for  the  rather  large  drag  C = 0.01  . 

The  flutter  speeds  fall  just  below  the  C = 0.01  divergence  speed. 

The  stability  boundaries  in  Fig.  6-13  do  not  allow  for  structional 
damping  and  give  no  feel  for  the  degree  of  stabilitv  at  speeds  near 
flutter.  To  gain  a better  understanding  of  the  type  of  instability 
that  has  been  found  with  steady  deformations,  and  also  define  the  over- 
all aeroelastic  behavior,  the  Laplace  transform  approach  detailed  in 
Chapter  V is  applied  to  the  large-  and  moderate-aspect-ratio  examples. 

The  true  aeroelastic  modes  are  then  conveniently  traced  in  the  complex 
plane  at  any  flight  conditions  using  root  locus  diagrams. 

Before  showing  the  root-locus  results,  it  is  interesting  to  see 
how  stability  is  suggested  bv  the  simple  harmonic  method  with  nonzero 
structural  damping  assumed.  The  C = 0 stability  boundary  of  Fig.  6-13 
is  reproduced  in  Fig.  6-18,  to  which  stabilitv  boundaries  for  three 
values  of  structural  damping  are  added.  The  sizeable  increase  in 
flutter  speed  with  g hints  that  the  instability  is  not  of  severe 
nature,  and  as  becomes  small  the  predominantly  fore-and-aft  bending 
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motion  is  obviously  very  lightly  damped.  For  small  deflections  with 
structural  damping  included,  there  is  a change  in  the  flutter  mode  back 
to  the  basic  bending-torsion  type  encountered  for  zero  si. ady  lift.  In 
spite  of  structural  damping's  stabilizing  influence,  there  still  exists 
the  possibility  of  a reduction  in  flutter  speed  at  high  load  factor 
flight  conditions. 

Root  locus  diagrams  depicting  the  dynamic  stability  of  the  moderate- 
and  large-aspect-ratio-wing  examples  with  steady  deflections  appear  in 
Figs.  6-19  and  6-20  respectively.  Dashed  lines  in  these  figures  show 
the  zero  steady  lift  loci  of  roots  for  increasing  speed,  which  are  the 
same  as  the  zero  drag  diagrams  of  Figs.  6-4 (c)  and  6-4 (g).  Solid  lines 
trace  the  elastic  modes  for  selected  constant  speeds  as  the  angle  of 
attack  a is  varied,  and  originate  for  a = 0 at  nonlifting  roots 
corresponding  to  these  speeds. 

In  addition  to  the  normal  modes  of  free  vibration  involving  vertical 

bending  and  torsion,  previously  seen  on  the  iaj  axis  in  Figs.  6-4,  the 

first  fore-and-aft  bending  mode  natural  frequency  now  also  must  be 

included.  This  normal  mode  remains  an  uncoupled,  undamped  aeroelastic 

mode  at  all  speeds  at  a = 0 , but  a family  of  constant  speed  branches 

emanates  from  this  root  with  steady  lift  included.  Since  its  natural 

frequency  u>  is  a factor  /F  larger  than  the  first  vertical  bending 
V1 

assumed  mode  frequency  to  , this  new  normal  mode  lies  on  the  iu 

W1 

axis  approximately  at  a multiple  /r  of  the  lowest  normal  mode  frequency. 
This  falls  above  the  zero-lift  flutter  frequency  for  the  moderate 
aspect  ratio  wing,  but  less  than  it  in  the  high-aspect-ratio  example. 
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A general  feature  of  Fig.  6-19  and  6-20  is  that  as  steady  lift 
increases  from  zero,  one  group  of  constant-speed  branches  tends  to 
stabilize  and  rapidly  increase  in  frequency,  whereas  another  lower 
frequency  family  of  root  paths  reduces  in  stability  and  decreases  in 
frequency.  This  latter  group  is  responsible  for  the  stability  boundaries 
of  Figs.  6-7  and  6-13.  In  the  high-aspect-ratio  example,  they  originate 
from  the  lowest  chordwise  bending  normal  mode,  but  for  the  moderate- 
aspect-ratio  wing  they  originate  on  the  zero-lift  torsion  branch  of  the 
root  locus.  The  constant-speed  branches  which  originate  on  the  first 
and  second  vertical  bending  zero  lift  paths,  for  both  examples,  do  not 
show  much  sensitivity  to  steady  deflections. 

To  illustrate  more  clearly  the  role  played  by  the  chordwise  bending 
upon  stability,  the  large-aspect-ratio  example  is  modified  by  increasing 
the  bending-stiffness  ratio  t from  25  to  60.  This  raises  the  dimension- 
less natural  frequency  of  the  first  chordwise  bending  mode  from  0.6216 
to  0.9629,  which  is  greater  than  the  zero  lift  bending-torsion  flutter 
frequency  of  0.8921.  The  stability  boundary  calculated  for  this  modified 
example  appears  in  Fig.  6-21  and  the  associated  flutter  modes  are 
presented  in  Figs.  6-22  and  6-23.  The  root  locus  obtained  via  the 
Laplace  transform  method  appears  in  Fig.  6-24(a)  and  the  true  stability 
of  constant-speed  branches  yielding  instability  is  better  depicted  in 
6-24 (b)  using  the  damping  ratio  C • 

The  discontinuity  in  the  stability  boundary  is  only  a consequence 
of  the  solution  procedure  of  Chapter  IV  and  covers  a region  where  solu- 
tions that  do  exist  could  not  be  determined.  It  is  due  to  an  interaction 
of  the  predominantly  second  bending,  stable  aeroelastic  mode  with  the 
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loots  that  yield  instability.  The  simple  harmonic  solution  method 
involves  fixing  the  angle  of  attack  ot  and  searching  for  the  neutrally 

**  .W  • «• 

stable  lifting  condition  whose  speed  matches  the  calculated  flutter 
speed.  In  Fig.  6-21,  the  gap  in  the  curve  falls  between  points  found 
for  a = .0095  and  .01  rad. 

For  all  points  on  the  left  segment  of  the  stability  boundary,  the 
simple  harmonic  solution  yielded  an  additional  highly  damped  eigenvalue 
whose  frequency  was  slightly  below  the  flutter  frequency.  On  the  right 
portion  its  frequency  was  above  the  flutter  frequency.  From  the  root- 
locus  diagram  (Fig.  6-24 (a))  this  highly  damped  eigenvalue  can  be 
identified  as  the  predominantly  second  vertical  bending  aeroelastic 
mode,  and  the  discontinuous  behavior  in  Fig.  6-21  coincides  with  the 
crossing  of  frequencies  as  the  downward  moving  constant  speed  loci 
associated  with  flutter  pass  the  stable  second  bending  aeroelastic  mode 
frequency.  By  correlation  with  the  sec  ion  VI-A  discussion  of  the  effect 
of  P and  C on  flutter,  the  general  effect  of  the  second  bending  mode 
appears  to  be  destabilizing  when  its  frequency  is  just  below  that  of 
flutter,  and  stabilizing  when  its  frequency  is  just  above. 

The  flutter  mode  shapes  (Fig.  6-22,  6-23)  reflect  an  interesting 
transition  as  flutter  frequency  drops  below  the  second  bending  aero- 
elastic mode  frequency  for  increasing  w . To  the  left  of  the 

o 

discontinuity  the  mode  shapes  resemble  the  zero-lift  large-aspec* -ratio 
flutter  behavior,  to  which  they  converge  as  a -*•  0 . For  large'  w , 
though,  the  modes  closely  resemble  flutter  mode  shapes  for  the  T = 25 
case  (Figs.  6-14,  6-15)  which  in  turn  resemble  mode  shapes  for  the 
moderate-aspect-ratio  example.  This  similarity  is  most  evident  in  q 

W2 

and  q . 
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In  Fig.  6-25  the  flutter  mode  shapes  are  diagrammed  for  two  lifting 
conditions  to  assist  physical  visualization  of  the  motion.  The  larger 
steady  deflection  illustration  is  similar  to  those  given  in  Fig.  6-16. 

But  the  smaller  steady  deformation  mode  shape,  with  its  large  contribu- 
tion by  the  second  assumed  bending  mode,  has  the  wingtip  rotating  in  a 
clockwise  direction  about  its  elliptical  path  for  a counterclockwise 
motion  at  midsemispan.  In  spite  of  the  relatively  small  steady  deflec- 
tion (a  deflection  of  1.226  semichords  at  the  wingtip  for  a semispan  of 
roughly  19  semichords)  the  motion  is  quite  three-dimensional,  indicating 
a significant  inertial  contribution  to  flutter  in  fore-and-aft  bending. 

The  three  root-locus  diagrams  (Figs.  6-19,  6-20,  and  6-24)  exhibit 
the  basic  effect  which  steady  deflections  have  upon  dynamic  stability 
when  incompressible  strip-theory  airloads  are  used.  The  basic  zero- 
steady-lift  bending-torsion  flutter  root  together  with  the  first  fore- 
and-aft  assumed  mode  produce  a pair  of  constant-speed  branches,  one  of 
which  rapidly  stabilizes  and  increases  in  frequency  while  the  other 
decreases  in  frequency  and  becomes  unstable  for  speeds  below  the  a = 0 
flutter  speed.  These  latter  aeroelastic  modes  are  lightly  damped,  and 
generally  the  onset  of  flutter  at  constant  speed  for  increasing  steady 
deflections  would  not  be  as  severe  as  for  that  encountered  with  increasing 
speed.  The  reduction  of  flutter  speed  with  steady  deformations  has  been 
noted  to  be  greatest  when  the  first  fore-and-aft  natural  frequency  is 
near  to  the  basic  bending-torsion  flutter  frequency. 

A better  comprehension  of  the  various  aeroelastic  modes  comes  from 
inspecting  mode  shapes  at  both  subcritical  and  supercritical  conditions. 

In  Fig.  6-26,  mode  shapes  for  selected  points  along  the  U = 6 and  7 
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branches  in  Fig.  6-19,  the  moderate-aspect-ratio  example,  are  shown  in 
phasor  form.  Likewise  for  the  high-aspect-ratio  wing  having  T = 60  , 
Fig.  6-27  presents  mode  shapes  for  certain  steady  deflections  along  the 
U = 6 and  7 paths  from  Fig.  6-24.  These  mode  shapes  disclose  that 
the  branches  which  stabilize  and  increase  in  frequency  consist  largely 
of  motion  in  the  first  assumed  torsion  mode,  having  a small  q 


1 


contribution  locked  in  a characteristic  phase  relationship  with  q of 


1 


near  180°.  The  branches  producing  instabilities,  furthermore,  have  q 


1 


nearly  in  phase  with  q and  give  the  previously  shown  flutter  mode 


1 


shapes  as  they  cross  the  iu>  axis.  Mode  shapes  (Fig.  6-27)  for  the 
essentially  second  bending  aeroelastic  modes  show  the  dominance  of  the 
second  assumed  vertical  bending  mode  in  this  branch,  which  remain  nearly 
fixed  in  the  complex  plane  at  a frequency  close  to  the  second  assumed 
bending  mode  natural  frequency. 

The  role  of  the  bending-stiffness  ratio  T in  flutter  of  lifting 
wings  is  next  examined.  Figures  6-28  and  6-29  give  flutter  speeds  and 
frequencies  found  for  the  same  moderate-aspect-ratio  example  used  earlier 
compared  with  solutions  for  different  values  of  x . For  T = 12  the 
dimensionless  natural  frequency  of  the  first  chordwise  bending  mode  is 
0.86124,  which  is  only  slightly  higher  than  the  zero-lift  flutter 
frequency  of  0.8332,  and  its  stability  boundary  shows  the  most  marked 


decrease  in  flutter  speed  as  w increases  from  zero.  The  other 

o 


extreme,  T = 1000  , has  a dimensionless  natural  frequency  in  chordwise 
bending  of  7.862,  vet  some  decrease  in  flutter  speed  still  occurs. 

Elastic  deformation  in  bending  about  the  major  principal  axis  of  the 
airfoil  cross  section  should  be  virtually  suppressed  and  all  fore-and-aft 
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motion  connected  with  vertical  bending-torsion  elastic  coupling.  Calcu- 
lations are  impractical  for  larger  T since  numerical  problems  begin 
to  appear,  such  as  convergence  difficulties  with  the  nonlinear  steady 
deflection  iterative  solution. 

The  uniform  decrease  in  flutter  frequencies  as  t is  decreased 

reflects  the  additional  inertia  from  the  larger  fore-and-af t motions 

brought  about  by  the  reduced  elastic  stiffness  in  bending  about  the 

airfoil  major  principal  axis.  Decreasing  flutter  frequencies  for 

increasing  steady  deformations  likewise  should  be  due  in  part  to  the 

increasing  contribution  of  q , although  the  increase  in  the  relative 

V1 

participation  of  q is  also  a factor. 

W1 

In  Fig.  6-30  the  stability  boundaries  shown  earlier  for  the  large- 
aspect-ratio  example  with  T = 25  and  60  are  compared  along  with 
curves  having  T = 10,  40,  and  200.  The  appearance  is  complicated  by 
the  interaction  of  the  flutter  modes  with  the  predominantly  second 
vertical  bending  aeroelastic  mode  as  discussed  earlier,  and  by  fore-and- 
aft  bending  natural  frequencies  sufficiently  low  to  cause  convergence 
as  a -*•  0 to  speeds  below  the  true  nonlifting  flutter  speed. 

For  both  T = 200  and  60  the  flutter  frequencies  descend  through 

the  range  of  the  second  vertical  bending  mode  frequency  as  w is 

o 

increased,  causing  the  discontinuous  stabilitv  boundaries.  The  flutter 
frequencies  (Fig.  6-31)  reflect  the  role  of  this  second  bending  mode, 
whose  dimensionless  natural  frequency  is  0.77904  and  (as  already  shown) 
a stable  aeroelastic  mode  with  approximately  this  frequency  exists  at 
flutter.  As  mentioned  in  discussing  the  T = 60  results,  these  curves 
are  actually  continuous,  but  the  simple  harmonic  solution  procedure 
could  not  produce  results  within  the  gaps. 
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For  the  x = 40  , 25  , and  10  cases,  stability  boundaries  converge 


j 


to  free  vibration  in  the  first  fore-and-aft  mode  as  a -*■  0 ; their 
frequency  curves  converge  to  the  respective  natural  frequencies  of  this 
mode.  In  the  case  T = 40  , this  frequency  is  0.7862,  quite  close  to 
the  second  vertical  bending  assumed  mode  natural  frequency  0.77904. 

This  near  resonance  caused  difficulty  in  precisely  locating  the  neutral 
stability  curves  for  small  w^  , as  the  damping  became  extremely  small; 
of  course  this  trouble  is  unimportant  since  even  small  structural  damping 
would  raise  flutter  speeds  considerably  here. 

The  stiffness  in  chordwise  bending  is  thus  a factor  in  the  flutter 
of  lifting  cantilever  wings.  Although  stability  is  most  adversely 
affected  for  near  the  zero-lift  flutter  frequency,  this  effect 

appears  over  a wide  range  of  I . Higher-aspect-ratio-wings  apparently 
experience  a greater  decrease  in  flutter  speed  with  steady  deformations, 

given  that  (i)  is  sufficiently  near  the  a = 0 flutter  frequency. 

V1 

The  mass  ratio  M is  the  only  parameter  which  can  change  for  a 
specific  wing,  as  it  depends  upon  altitude.  In  Fig.  6-32  the  dependence 
of  both  flutter  and  divergence  speeds  on  M for  the  moderate-aspect- 


ratio  example  appears  for  the  complete  range  of  mass  ratios  of  practical 
interest.  Divergence  speeds  for  several  values  of  steady  drag  together 
with  zero-lift  flutter  speeds  for  C = 0 and  0.01  are  compared.  In 
addition  zero  drag  flutter  speeds  for  steady  lift  giving 


are  shown  for  M = 10  , 40  , and  100  . This  same  information  is 


depicted  in  Fig.  b-33  for  the  large-aspect-rat io  example  except  wfngtip 
deflections  of  2 and  4 semichords  are  used. 

Divergence  speeds  for  all  drag  values  vary  exactly  as  rrl  , as  can 
be  seen  from  the  divergence  determinant  (2-29),  where  M appears  only 
as  a product  with  the  inverse  square  of  divergence  speed.  The  0 = 0 

| ■ 

divergence  speeds  for  these  two  examples  are  identical,  but  11  decreases 
more  rapidly  with  C for  the  higher  aspect  ratio. 

For  mass  ratios  above  about  3,  nonlifting  flutter  speeds  behave 
approximately  as  . A minimum  in  flutter  speed  is  found  near  M 3 . 

whith  a rapid  asymptotic  rise  to  infinity  following  a further  decrease 
in  M . This  reflects  well  known  results  for  incompressible  flow  (Ref.  3, 
page  247).  A practical  application  for  mass  ratios  sufficiently  small 
to  be  theoretically  free  from  flutter  is  the  stability  of  hydrofoils 
used  in  high-speed  marine  transportation.  Divergence  would  be  the  type 
of  instability  encountered.  This  conclusion  emphasizes  the  importance 
of  allowing  for  steady  drag  since  its  effect  on  divergence  is  propor- 
tionally the  same  for  any  mass  ratio  and  the  decrease  of  divergence  speed 
with  C can  be  considerable. 

Divergence  is  more  important  relative  to  flutter  for  mass  ratios 
around  10  than  for  higher  mass  ratios;  this  indicates  a greater  likelihood 
of  light  aircraft  and  sailplanes  experiencing  divergence.  Calculation 
of  divergence  speeds  for  a high-performance  sailplane  appears  as  a likely 
application  in  which  allowance  for  steady  drag  effects  would  be  desireable. 

The  effect  of  steady  deformation  on  flutter  appears  insensit ive  to 
M . For  equal  tip  deflections,  0 reduces  hv  about  the  same  proportion 
at  each  of  the  three  mass  ratios  in  the  t wo  figures.  For  the  high-aspect- 
ratio  wing,  these  deflected  flutter  speeds  correspond  I o the  aeroelastic 
mode  connected  with  chordwise  bending,  as  discussed  for  M 

i r\ 
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The  effect  of  changes  tn  the  remaining  parameters  A , S , and  i 
has  also  been  investigated  for  the  steady  deformed  case,  but  no  new 
interesting  behavior  was  found.  Their  influence  at  zero  lift,  as 

•i 

tabulated  in  Ref.  1,  is  carried  over  into  the  deflected  flutter  behavior. 
Stability  boundaries  and  flutter  mode  shapes  for  increasing  steady  lift 
show  the  same  alterations  of  the  basic  flutter  phenomenon  as  is  demon- 
strated by  the  examples  used  here. 

To  clarify  the  strong  influence  which  steady  chordwise  loads  can 
have  upon  divergence,  particularly  for  high  aspect  ratios,  the  dependence 
of  Up  upon  C as  aspect  ratio  parameter  P is  varied  appears  in 
Fig.  6-34.  The  corresponding  zero-lift  flutter  speeds  for  this 
configuration  (Fig.  6-1)  can  be  compared  directly.  Unquestionably,  for 
higher  aspect  ratios  (P  < 0.01)  drag  forces  typically  encountered  in 
flight  (i.e.,  C = 0.0023)  can  cause  divergence  to  become  no  less 
important  for  flight  safety  than  primary  bending-torsion  flutter. 

C . Effect  of  Unsteady  l.e  ading-Edgo  Suction  Forces  from  Two  Pinions  io  tu  i J 
Inc o mpre ss i hie  Flow 

Substituting  the  terms  in  (5-20)  for  their  counterparts  in  the  aero- 
dynamic matrix  (5-3)  allows  inclusion  of  the  linearized  unsteady  chord- 
wise  forces  arising  from  incompressible  potential  strip  theory  as 
described  in  Chapter  V.  This  effect  is  present  only  for  a t 0 and 
will  grow  as  angle  of  attack  is  increased.  The  initial  trend  of  stability 
as  oi  increases  from  zero  should  be  unchanged  from  the  results  already 
discussed,  but  at  higher  steady  lifting  conditions  leading  edge  suction 
should  have  a noticeable  influence. 
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Figure  6-35  displays  how  allowance  for  the  linearized  unsteady 
propulsive  force  affects  the  destabilizing,  constant  speed,  increasing 
steady  lift  branches  of  the  root  locus  in  Fig.  6-19,  which  is  the 
moderate-aspect-ratio-wing  example.  For  this  case  a tip  deflection  of 
just  one  semichord  is  a very  large  deformation  assuming  a conventional 
bending-to-torsion  stiffness  ratio  EI^/GI^  , since  the  semispan  would 
be  about  9 semichords.  Thus  larger  deflections  shown  are  only  of 
academic  interest,  since  they  exceed  in  practice  the  limits  imposed  by 
several  underlying  assumptions. 

Figure  6-36  recasts  the  Fig.  6-35  information  into  a better  form 
for  inferring  stability;  the  damping  ratio  of  the  same  root  branches  is 
plotted  for  wingtip  deflection.  Stability  is  slightly  increased  when 
leading-edge  suction  is  added,  but  for  practical  deformations  its  effect 
is  really  not  too  significant.  A stability  characteristic  of  deformed 
wings,  also  evident  in  the  format  of  Fig.  6-36,  is  the  reduced  severity 
of  flutter  onset  with  increasing  speed  at  constant  wQ  for  higher  load 
factors.  | 

Figure  6-37  gives  the  manner  in  which  the  root  locus  for  the  high- 
aspect-ratio  example  (Figure  6-20)  is  altered  by  allowing  for  linearized 
unsteady  suction  forces.  Shown  is  the  destabilizing  family  of  constant- 
speed  branches  associated  with  the  first  fore-and-aft  bending  normal 
mode.  Again  true  stability  with  steady  wingtip  deflection  is  also 
shown  (Fig.  6-38).  Also  added  are  the  constant-speed  branches  which 
stabilize  and  increase  in  frequency  as  steady  deformations  increase. 

Once  again  the  unsteady  suction  force  has  a slightly  stabilizing 
influence.  Given  the  same  ratio  EI^/GI^  , this  example  has  twice  the 
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span  in  semichords  of  that  for  the  moderate-aspect-ratio  example,  hence 
larger  wingtip  deflections  in  semichords  can  be  tolerated.  After 
allowing  for  this  factor,  the  unsteady  chordwise  potential  forces  do 
not  appear  to  have  any  greater  influence  on  stability  for  the  larger 
aspect  ratio. 

Figure  6-39  shows  how  stability  of  the  modified  high-aspect-ratio 
example  (having  T = 60  , presented  in  Fig.  6-24)  depends  on  steady 
deformations  and  on  the  linearized  unsteady  propulsive  force.  Comparison 
with  Fig.  6-38  shows  a similar  effect  due  to  leading-edge  suction,  but 
steady  deformations  definitely  have  a stronger  and  more  immediate 
destabilizing  effect  for  the  stiffer  chordwise  bending  case.  This  is 
due  to  the  proximity  of  the  fore-and-aft  bending  normal  mode  frequency 
to  the  zero-lift  bending-torsion  flutter  frequency. 

Inclusion  of  the  linearized  unsteady  leading-edge  suction  terms 
derived  in  Section  D of  Chapter  V is  thus  found  to  be  stabilizing  for 
deformed  wings.  Aeroelastic  modes  involving  a primary  contribution  by 
the  first  chordwise  bending  mode  appear  to  be  most  affected  by  these 
terms.  With  their  inclusion,  potential  flow  strip  theory  has  been 
fully  exploited  for  this  problem,  and  further  improvement  in  the  aero- 
dynamics involves  the  compressible  three-dimensional  loads  of  Chapter 
VII. 


D . Two  P ractlc a 1 F.xamp  1 es 

The  manner  in  which  chordwise  forces  and  steady  deformations 
influence  the  aeroelastic  stability  of  cantilever  wings  has  thus  lar 
been  illustrated  by  means  of  idealized  examples.  After  one  has  identi- 
fied the  fundamental  effects,  it  is  instructive  to  apply  the  same 
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solution  techniques  to  typical  designs.  Sailplanes,  having  large 
aspect  ratios,  low  mass  ratio,  and  low  operating  speeds  (which  permit 
the  use  of  incompressible  aerodynamics),  are  a logical  choice  for 
practical  examples.  Accordingly,  two  existing  sailplanes  are  modeled 
within  the  approximations  imposed  by  the  assumption  of  uniform  spanwise 
mass  and  stiffness  properties  and  strip  theory  airloads.  Their  stability 
is  studied  using  the  same  techniques  applied  before. 

The  first  example  is  modeled  from  information  given  about  the 

Slingsby  Dart  17R  in  Ref.  15,  as  summarized  in  Table  6.1.  The  approach 

taken  in  determining  the  wing  stiffnesses,  in  lieu  of  looking  at  wing 

construction  details,  was  to  use  photographs  in  Ref.  15,  from  which  tip 

deflection  and  twist  at  a load  factor  of  4 could  be  measured.  The 

stiffness  EI^  so  obtained  is  clearly  larger  than  would  be  expected 

for  the  weight  and  type  of  construction;  this  inaccuracy  is  due  to  the 

uniform  stiffness  restriction,  the  rectangular  planform,  and  strip  theory 

steady  airloads.  The  torsion  stiffness  GI,  is  more  reasonable.  The 

a 

ratio 


turns  out  impractically  large.  This  number  leads  to  a larger  value  of 

P than  would  be  expected,  which  causes  bending  natural  frequencies  to 

be  larger  relative  to  torsion  natural  frequencies  than  would  probably 

occur  on  the  actual  vehicle.  It  can  be  added  that,  had  GI,  been 

d 

increased  to  compensate  for  the  seemingly  excessive  EI^  , then  the 
flutter  speeds  would  have  been  absurdly  high.  Since  wing  construction 
details  were  not  available,  typical  values  for  the  parameters  A , S, 
i^  , and  T were  simply  assigned. 
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A different  philosophy  was  applied  when  modelling  the  second 
vehicle,  the  Gemini  two-place  high  performance  sailplane  detailed  in 
Ref.  16.  With  construction  details  of  the  metal  wing  structure  described 
in  Ref.  lb,  it  was  possible  to  assign  reasonably  accurate  values  to  all 
six  of  the  dimensionless  parameters  from  a diagram  of  the  average  wing 
cross  section  (Fig.  6-40).  The  pertinent  details  of  modeling  the 
Gemini  are  listed  in  Table  6-2. 

One  notable  result  is  a much  lower  bending  stiffness  EI^  than 
for  the  Dart  17R.  This  result  will  compromise  the  static  deflection- 
load factor  relationship  but  should  favor  the  dynamic  modeling.  The 
considerable  difference  in  bending  stiffnesses  between  the  two  examples 
is  revealed  in  Fig.  6-41,  which  shows  how  the  true  load  factor  (found 
using  tlie  respective  aircraft  gross  weights)  varies  with  steady  vertical 
tip  deflection  at  speeds  near  the  expected  flutter  speeds.  The  poor 

|1 

static  modeling  of  the  Gemini  should  not  adversely  affect  its  flutter 
results,  however,  tip  deflection  rather  than  load  factor  should  be  used 

* 

to  measure  the  amount  of  steady  lift. 

The  mass  per  unit  span  of  the  Gemini  wing  was  taken  to  be  the 
average  for  the  outer  two-thirds  of  the  semispan,  so  as  better  to  model 
it  dynamically.  The  mass  ratio — over  twice  that  for  the  Dart  1 7R — 
reflects  the  heavier  construction  needed  for  a two-place  vehicle  and 
the  smaller  average  semichord.  This  design  (Ref.  16)  intentionally 
has  a higher  wing  loading  to  optimze  the  glide  slope,  with  thermalling 
performance  Improved  using  full-span  flaps.  The  chordwise  bending 
stiffness  was  conservatively  estimated  yet  still  yielded  a quite  low 
value  of  r = 5 . The  low  EI^/GI^  ratio,  together  with  the  very  large 


138 


aspect  ratio,  result  in  the  extremely  small  P = 0.00128  . This  wing 
thus  is  an  extreme  case  in  the  context  of  the  examples  studied  earlier. 

The  effect  of  steady  drag  on  divergence  speeds  and  nonlifting 
flutter  speeds  is  shown  in  Fig.  6-42  for  the  Dart  17R  example  and  in 
Fig.  6-43  for  the  Gemini.  For  both,  the  divergence  speed  drops  below 
the  zero  lift  flutter  speed  at  drag  values  which  can  be  realistically 
expected  in  flight.  The  excessively  large  P of  the  Dart  17R  model 
moderates  the  drop  in  its  divergence  speed  due  to  drag.  For  flutter 
only  a slight  dependence  of  speed  on  C is  found.  On  the  Gemini,  in 
fact,  the  C = 0.02  flutter  speed  is  only  0.28%  less  that  for  C = 0 . 
Its  larger  S tends  to  lower  the  flutter  speed  relative  to  divergence. 

Figure  6-44  depicts  the  stability  boundary  as  affected  by  steady 
wingtip  deflection  for  the  Dart  17R.  The  associated  flutter  mode  shape 
amplitudes  and  phase  relations  appear  in  Figs.  6-45 (a),  (b).  The 
coupling  effect  between  the  first  fore-and-aft  bending  mode  and  the 
primary  zero-lift  bending-torsion  flutter  mechanism  appears  to  be 
responsible  for  a decrease  in  flutter  speed  with  lift,  as  in  earlier 
examples.  The  corresponding  root  locus  diagram  (Fig.  6-46)  reveals  that 
the  first  chordwise  bending  natural  frequency  by  coincidence  happened 
to  fall  almost  exactly  on  the  zero  lift  flutter  frequency.  This 
phenomenon  causes  the  pronounced  drop  in  flutter  speed  for  small  deflec- 
tions seen  in  Fig.  6-44.  Inclusion  of  the  unsteady  propulsive  force  is 
slightly  stabilizing  at  large  deflections,  as  in  previous  cases. 

A parallel  analysis  of  the  Gemini  example,  given  by  Figs.  6-47, 
6-48,  and  6-49,  yields  a very  different  response  to  steady  deflections. 
Flutter  speed  decreases  only  slightly  at  representative  tip  displace- 
ments, and  the  mode  shape  contributions  from  q , q , and  q,  show 


r- 

almost  no  sensitivity  to  w^  . The  chordwise  bending  contribution 
suggests  the  cause,  made  clear  in  the  root  locus  diagram.  The  combina- 
tion of  a very  small  P and  small  t produces  a fundamental  fore-and-aft 
bending  natural  frequency  which  is  much  smaller  relative  to  the  zero 
lift  flutter  frequency  than  in  any  example  yet  treated.  Furthermore, 
the  second  chordwise  frequency  drops  to  near  the  flutter  frequency. 

This  second  mode  (rather  than  the  fundamental)  participates  in  flutter 
of  the  deformed  wing,  yet  it  is  not  as  strongly  coupled  elastically  with 
the  first  vertical  bending  and  torsion  modes.  The  aeroelastic-mode 
branches  associated  with  the  first  chordwise  bending  mode  no  longer 
play  an  important  role  in  stuuility  with  steady  deformations.  Indeed 
it  has  degenerated  into  a virtually  uncoupled,  neutrally  damped  fore-and- 
aft  vibration.  The  linearized  effect  of  unsteady  leading-edge  suction 
on  this  type  of  flutter  is  negligible. 

Evidently,  for  extreme  cases  such  as  the  Gemini  model  with  a fore- 
and-aft  bending  frequency  much  lower  than  the  zero-lift  bending-torsion 
flutter  frequency,  the  lowest  chordwise  mode  will  not  participate  strongly 
in  flutter.  The  low-frequency  root  branches  associated  with  it  give  no 
cause  for  concern.  Before  this  calculation  is  accepted  as  a definitive 
flutter  analysis  of  an  existing  sailplane,  however,  the  cantilever  root 
boundary  conditions  must  especially  be  recalled.  While  these  results 
do  offer  insight  into  symmetric  flutter  of  the  actual  vehicle,  the 
possibility  of  anti-symmetric  motions  involving  fuselage  roll  and  yaw 
is  entirely  suppressed.  For  the  Gemini  in  particular,  there  remains 
the  likelihood  that  anti-symmetric  vertical  and  fore-and-aft  bending 
modes  may  couple  to  produce  a different  type  of  flutter — possibly  one 
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with  a greater  sensitivity  to  steady  deformations.  Although  the  results 
here  are  insensitive  to  both  steady  deflections  and  to  chordwise  forces, 
steady  drag  unquestionably  plays  a critical  part  in  divergence. 


DART  17R  SPECIFICATIONS  (REF.  15): 


SPAN 

AREA 

ASPECT  RATIO 

GROSS  WEIGHT 


55.8  ft. 
149  ft2 

20.9 
800  lb. 


ASSUMPTIONS: 


WING  WEIGHT 
AVG.  SEMICHORD 


200  lb  =>  m - 0.1242 
1.33  ft. 


FROM  REF.  15  PHOTOGRAPHS: 

(wo)tip  =■  3 ft.  at  4g  =>  EIx  ~ 1,500,000  lb  ft' 
($0)TIp  = 1-5°  at  4g  *>  GId  = 340,000  lb  ft2 

FROM  THIS  INFORMATION  ONE  CAN  SPECIFY 

P ~ 0.01  M ~ 9.4  (sea  level) 

ASSUME  THE  REMAINING  PARAMETERS 

i = 0.25  S - 0.1 

a 

A “0.1  t = 25. 

EXPRESS  SPEED  V IN  FT. /SEC.  IN  TERMS  OF  U : 


l /J_ 

b /GI. 


V = 0.008432  V 


LOAD  FACTOR  AT  SEA  LEVEL: 


4HpV2biL  . ? 2 1 

" ■ (“  + n UFIT  q*1) 


TABLE  6.1  Modeling  of  the  Dart  17R  Wing 


GEMINI  SPECIFICATIONS  (REF.  16) 


ASPECT  RATIO 

GROSS  WEIGHT 

TOTAL  WING  WEIGHT 

WEIGHT  OF  OUTER 
20  FT. 


60.5  ft. 
124  ft2 


1065  lb. 
400  lb. 


110  lb.  =>  m = 0.1708  -^S- 

ft. 


ASSUME  b = 1.00  ft. 

PROPERTIES  ESTIMATED  FROM  THE  TYPICAL  SECTION  (Fig.  6-40) 
EIx  444,900  lb-ft2 

GId  402,400  lb-ft2 

El  2,179,000  lb-ft2 

J 0.04697  slug-ft. 


0.2322  ft. 


RESULTING  DIMENSIONLESS  PARAMETERS 


M =22.9 

P = 0.00128 

i = 0.275 

a 


A = 0.1 
S = 0.23 
T = 5 


EXPRESS  SPEED  V IN  FT. /SEC.  IN  TERMS  OF  U 


u = r /-FT-  v = 0.010335  V 
D v GI  , 
a 
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Re(p) 


FIGURE  6-4 (f) 
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BENDING 

BRANCH 


U = 4. 


FIGURE  6- 


5 Mode  Shapes  for  Selected  Roots  on  Each  Branch  of  Fig.  6-4 (g) 
for  C - 0 


FIGURE  6-6(a)  Flutter  S 
C as  P 
Mass  Ratli 


AMPLITUDE 


360. 


.001  .01  .1 


P 

FIGURE  6-6 (b) 
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Flutter  Speeds  as  Affected  By  Steady  Deformations;  Large-Aspect-Ratio  Exampl 


PHASE  DEG. 


r until 


-is 


Flutlfr  Mode  Shupt  nh,wi'  Irltl  lorn  fur 
(nrrr«t>«fldlnr  !*•  ff (.  *•  II 


Vrn  of 


'’a 


1 


path  of  wingtip 

path  of  section  at  y =.5 

FIGURE  6-16  Physical  Appearance  of  Flutter  Mode  of  Deflected  High 
Aspect  Ratio  Wing  of  Fig.  6-13 


FIGURE  6-17  Comparison  of  Drag  Effect  Upon  Divergence  and  Flutter  At  Zero  Steady  Lift  for  the 
Moderate-  and  High-Aspect-Ratio  Examples 


6-13  Stability  Bounda 


UNCLASSIFIED  SUDAAR-508  AFOSR-TR-7S-0856  NL 


3 of  3 


AO 

AQ53640 


constant  U, 
increasing  w 

zero  w„ , 
increasing  U 


-.4 

1 

W 

1 

Io 

• 

1 0.  + 

Re(p) 

Figure 

6-24(a)  Locus  of  Roots  for  True  A< 

;roelastic  Modes,  Large- 

Aspect-Ratio  Example  with 

T = 60 

.66  1.29 


FORE-AND-AFT 

q 1*,  BENDING 

BRANCH 

U = 6. 

2.44 


FIGURE  6-26  Mode  Shapes  for  Selected  Roots  from  Fig.  6-19 
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FORE-AND-AFT 

BENDING 


(flutter  mode) 


{ B°}tip  - 1.31  2.57  5.29 


FIRST 
VERTICAL 
BENDING 
BRANCH 
U=  6. 


FIGURE  6-27(a)  Mode  Shapes  for  Selected  Roots  from  Fip.  6-24 
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FIGURE  6-29  Effect  of 


On  Flutter  Speeds,  Large-Aspect-Ratio  Example 


FIGURE  6-31  Effect  of  T On  Flutter  Frequencies,  Large-Aspect-Ratio  Example 


40.  60.  80.  100. 

M 


Figure  6-32  Effect  of  Mass  Ratio  on  Divergence  and  Flutter,  Moderate- 
Aspect-Ratio  Example 


f 


85 


Re  ( p) 


FIGURE  6-35  Effect  of  Including  Unsteady  Leading  Edge  Suction  Forces 
Upon  the  Constant  Speed  Branches  of  Fig.  6-19 
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FIGURE  6-36  Damping  Ratio  of  the  Branches  in  Fig.  6-34 


LOAD  FACTOR 


6 


1 


FIGURE  6-41  Comparison  of  Vertical  Bending  Deflections  of  the  Two 
Sailplane  Wing  Models  for  filven  Steadv  Load  Factor 
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DIVERGENCE 

SPEEDS 


FIGURE  6-44  Stability  Boundary  For  the  Dart  17R  Model 
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I woItIP  (feet) 


FIGURE  6-48  Flutter  Mode  Shape  Amplitudes  and  Phase  Angles  for  the 
Gemini  Model,  for  Unit  Amplitude  and  Zero  Phase  of  q. 


Chapter  VII 


FLUTTER  VELOCITY  USING  AIRLOADS  FROM  THREE-DIMENSIONAL 
SUBSONIC  AERODYNAMIC  THEORY 

A.  Inclusion  of  Three-Dimensional  Airloads 

The  influence  of  steady  deformations  and  chordwise  forces  upon 
dynamic  stability  of  the  uniform  cantilever  wing  has  been  examined  in 
Chapter  VI  using  lifting  airloads  predicted  by  incompressible  steady 
and  unsteady  strip-theory.  This  ar proximate  modeling 

of  the  aerodynamic  loads  made  possible  their  convenient  numerical  compu- 
tation for  any  convergent,  neutral,  or  divergent  oscillations  of  interest. 
As  a result,  the  iterative  solution  schemes  of  Chapters  IV  and  V could 
be  developed  and  a variety  of  wing  configurations  could  be  analyzed 
efficiently. 

The  accuracy  with  which  the  incompressible  strip  theory  results 
approximate  the  three-dimensional  compressible  flow  situation  is  next 
explored  by  extension  to  subsonic  three-dimensional  lifting  airloads. 

The  simple  harmonic  flutter  solution  method  described  in  Chapter  IV  is 
accordingly  modified  to  use  subsonic  three-dimensional  steady  and 
oscillating  unsteady  aerodynamic  loads  calculated  separately  by  the 
computer  program  written  by  Rowe,  et  al.  (Ref.  18).  Results  are  then 
found  which  demonstrate  that  the  phenomena  discussed  in  Chapter  VI 
still  occur  after  three-dimensional  aerodynamics  are  introduced.  The 
role  of  unsteady  potential  chordwise  loads  in  flutter  is  investigated 
as  is  the  effect  of  compressibility.  Since  the  use  of  externally 
computed  air  loads  requires  a more  cumbersome  solution  method,  only 
enough  results  are  sought  to  provide  direct  comparison  with  the  incompres- 
sible strip-theory  calculations. 


i 
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The  Rowe  computer  program  solves  the  pressure-downwash  Integral 
equation  (Ref.  24)  for  compressible  flow  about  a steady  or  oscillating 
planform.  The  kinematic  downwash  boundary  condition  for  each  structural 
mode  is  enforced  by  collocation,  at  a set  of  user-specified  points,  of 
the  downwash  distributions  associated  with  an  assumed  series  of  pressure 
functions.  For  this  application  seven  collocation  chords  are  specified 
having  five  collocation  points  per  chord.  Six  elastic  structural  modes 
are  input — three  in  vertical  bending  and  three  in  torsion.  The  three 
torsion  modes  are  defined  for  A = .1  (lateral  displacements  depend  on 
elastic  axis  location).  A thorough  description  of  the  theoretical 
aspects  of  the  subsonic  kernel  function  program  is  provided  in  Ref.  18, 
and  the  programming  details  are  documented  in  Ref.  19.  Its  capability 
for  modeling  trailing-  and  leading-edge  control  surfaces  is  not  required 
in  the  present  application. 

The  unsteady  potential  chordwise  forces  can  be  deduced  from  program 
output,  as  follows.  For  each  structural  mode  and  its  downwash  the 
program  calculates  the  complex  amplitude  of  an  associated  distributed 
lifting  pressure  difference  on  the  rectangular  planform,  which  in 


dimensionless  form  for  the  jth  mode  is,  per  unit 
simple  harmonic  time  dependence. 


qj 


having  steady  or 


AP  (x,y) 

(7-1)  ACPj(x,y)  = j^pVa 


? z ,«>£<‘>(y)8('‘)(x) 
V*1  y=l  y 


Here  a^  are  coefficient  multipliers  of  the  series  expansion  of 
pressure  on  the  planform,  which  can  be  listed  in  the  program  output. 
The  assumed  spanwise  pressure  distributions  are 
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r 

< 


(7-2) 


f(v)(y) 


3ln(2v-l)$ 

ain't1 


v - 1,2,3, ...N 


(7-3)  <t>  - cos  1(^) 

N - number  of  downwash  chords  on  semispan 

The  assumed  chordwise  pressure  distributions,  dependent  upon  x only 
for  a rectangular  planform,  are 

(7-4)  g^(x)  - / cot  U - l 

I sin(u-l)0  U = 2,3,4, ...M 

(7-5)  0-cos_1(-J) 

M = number  of  downwash  points  on  a downwash  chord 

where  x is  measured  aft  from  the  midchord. 

The  resultant  chordwise  component  of  potential  airloads,  acting 
in  the  positive  x direction,  can  be  expressed  as 

, 3z 

(7-6)  D (y;t)  = - 'spV2  / AC  (x,v;t)  (x,v;t)dx  - F (v;t) 

P *“D  p C/  X S 

The  first  term  represents  the  x-component  of  the  force,  which  is  normal 
to  the  deflected  chord,  second  term  contains  the  contribution  of  leading- 
edge  suction.  It  is  an  idealization  of  linearized  theory,  which  is 
supposed  to  approximate  the  actual  effects  of  low  pressure  acting  around 
a curved  leading  edge. 

Steady  and  unsteady  parts  of  the  pressure  and  mean-surface  chordwise 
slope  can  be  separated: 


J 
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The  suction  force  likewise  can  be  separated.  The  leading-edge  inverse- 
square-root  pressure  singularity  strengths  for  steady  and  unsteady  flow 
can  be  defined, 

(7-10)  CF  (y)  = ^ lim  [-/x+b  AC  (x,y)] 
o x-*— b o 

(7-11)  C (y;t)  = 7 lim  [-  /x+b  AC  (x,y;t)] 

F1  4 ^-b  P1 


From  equations  (5-9)  and  (5-10),  the  leading  edge  suction  force  in 
terms  of  the  vorticity  singularity  is 

(7-12)  F (y;t)  = £ p[  lim  (v^+b  Y(x,y;t))]2 

a x+-b 

where  the  effect  of  compressibility  is  now  included  with  the  /i_m2 

a 

factor  (deduced  from  Eq.  12-1  of  Ref.  26).  Vorticity  and  pressure  dis- 
continuity distributions  in  the  vicinity  of  the  leading  edge  are  related 


1 


I 

I 


i 
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by  (Ref.  23,  Eq.  (5-93)) 


(7-13) 


ACp(x,y;t)  = 


- -y^(x,y;t) 


The  suction  force  in  terms  of  the  leading  edge  pressure  singularity 
is,  therefore. 


(7-14) 


F (y;t)  = /i-m2  np{X  Um  f-  ✓x+b  AC  (x,y;t)]}‘ 
s a 4 P 

x-1— b 


= A-M2  np{c  (y)  + C (y;t)}" 
a Fo  F1 


= /l^m2'  np[c^  + 2c  c + c;  ] , 

a Fo  Fo  F1  F1 


where  (7-7),  (7-10),  and  (7-11)  have  been  used.  Equations  (7-9)  and 
(7-14)  suggest  that  chordwise  forces  do  affect  both  the  steady  displace- 
ment solution  and  the  linearized  unsteady  stability  problem. 

Actual  computation  of  the  suction  force  contributions  can  be 
accomplished  through  program  output  of  the  series  coefficients  a^ ^ . 
Insertion  of  (7-1)  into  (7-10)  and  (7-11),  for  the  jth  mode,  leads  to 
the  steady  and  unsteady  leading-edge  singularity  strengths 


(7-15) 


:(J)(y)=-2v[^  ? a4<,3)f(v)(y)]^ 


* v-l  vl 


(7-16)  c£j)(y)  = - 2V[v^Z  z f(v)(y)]v^ 


Here,  in  taking  the  limit  x -*■  -b  , only  the  U = 1 chordwise  pressure 


distribution  terms  from  (7-4)  remain  since 


A 


(7-17)  lim  [ /x+b  g 
x-*~  b 


(U) 


(x)  i = r u = i 

L 0 H > 1 


The  single  summations  that  define  the  modal  spanwise  leading-edge 
pressure  singularity  strength  can  be  arranged  for  computational  purposes 

(7-18)  F°(y)  = v— ^ Z f(V)(y) 

3 v=l  o 

(7-19)  F (y)  = l a(j ) f(v>(y) 

1 * v=l  V11 


This  notation,  together  with  summation  over  all  modes,  allows  the  total 
singularity  strengths  (7-10)  and  (7-11)  to  be  computed  by 


(7-20)  Cp  (v)  = - 2V^b  E F*(y)q° 
Fo  j j j 


(7-21)  C (y , t)  = - 2V/b  E F (y)  q (t) 
1 .1  J -1 


Insertion  of  (7-20)  and  (7-21)  into  (7-14)  gives  for  the  suction  force 


(7-22)  F (y , t)  = 4IIpV2b  /wF  {[E  F°(y)q°r 

S V V 


+ 2 E E F* (y)q°  q,elu,tr  (y) 


v j 


V'"MV  J' 


■ i 


+ [E  F (y)  q elu)tl2} 
j 3 3 


The  notation  of  ( 7— J 8)  and  (7-19)  is  next  adjusted  to  identify  the 
specific  structural  modes  involved.  For  the  steady  problem  a rigid 
pitching  mode  must  be  used  to  solve  for  pressures  and  loads  due  to 
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airplane  angle  of  attack,  and  Its  associated  singularity  strength 


parameter  is  denoted  as  F^  . For  the  jth  elastic  torsion  mode  F® 

is  used;  of  course,  the  vertical  bending  modes  introduce  no  steadv 

lifting  loads.  Unsteady,  leading-edge-singularity  strength  parameters 

for  the  jth  elastic  torsion  and  bending  modes  are  specified  by  F^ 

and  F , respectively.  The  computational  form  of  (7-22)  is  then 

W1 

(7-23)  F (y;t)  = AIIpV2/!^  {[F°a+  E F"  q°  ]z 
s a a v=1  <pv  <J>V 


n n w , , 

+ 2 E [F#a  + E F®  q°  IF  -r1  e1Wt 
, , a , d>  M(})  w,  b 

j=l  V=1  % j 


n n 

r rr°~  j.  v t°  _ o 


+ 2 E [F  a * I f;  .•  IF  q,  e‘ 

j-1  “ v-1  t V *J  *j 


n w n - r\ 

+ f(E  F -r1  + E F q,  )e‘  ]2} 

J-1  “j  b .1-1  *1  *1 


The  x-component  of  the  ACp-force,  as  indicated  in  (7-9),  must  also 
be  expressed  in  a form  which  permits  computation.  Insertion  of  the 
modal  quantities 


(7-24)  *o  ■ £ f 

V * 1 


(7-23)  ACp  = ACp  a + E ACp.  q° 
o a , <p  m 

V=1  rv  V 


(7-26)  ^ ^ 


w.  n 


(7-27)  ACp  = £ ACp  -r1  + £ ACp.  q. 

1 j-i  ”j  b j-i  *j  *i 
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into  (7-9)  leads  to 


(7-28) 


D (y;t)  - *iPV2{[  E / ACp“  f dx  a qj 
V V-l  ^v 


+ ”,  V-b4Cp« f* dx ■>;  ■>; 1 - 1 ” " f*  *■  i;  -r 

U=1  v-1  vv  vv  i=l  v-1  i vv  vv 

+ E ^ /bbACp  f dx  qj  q 

1=1  v-l  i v ri 

n , n n . 

+ £ / ACp  f dx  a q + E E / ACp  f dx  q°  q,  ]e  Wt 

i-1  _b  a *i  *i  i-1  v-1  "b  *v  *i  X *i 


n n "w. 

+ [ E E / ACp  f , dx  -r—  q 
i-lj-l-b  wi  *j  b *] 

+ E E /b  ACp , f dx  q,  q ]e2iaJt  . 

i-1  j-1  "b  *1  *S  *i  *.1  - Fs(y;t) 


where  steady,  linear  unsteady,  and  nonlinear  unsteady  terms  in  the 
generalized  displacements  are  grouped  together. 

The  chordwise  integrals  in  (7-28)  involving  the  steady  and  unsteady 
modal  pressures  can  be  computed  at  any  spanwise  station  by  taking 
advantage  of  the  Rowe  program's  capability  to  compute  sectional  general- 
ized forces.  Program  output  is  of  the  form 

(7-29)  ^(y)  = -(lb  ACpj  (x,y)Hi(x,y)dx 

where  the  deflection  mode  shape  H^x.y)  for  all  integrals  in  (7-28) 

will  be  f.  (v)  . 

♦i 

The  nonlinear  modal  equations  (4-9)  for  steady  displacements  must 
be  modified  by  substitution  of  the  steady  compressible  lifting  airloads 
in  place  of  the  strip  theory  loads  (4-3).  Inclusion  of  the  steady 
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parts  of  (7-28)  and  (7-23)  will  account  for  the  induced  drag  caused  by 

the  three-dimensionality  of  the  flow.  Generalized  forces  which  will 

appear  in  equations  (4-9a)  and  (4-9c)  are  to  be  computed  directly,  but 

the  generalized  forces  for  the  chordwise  bending  equations  (4-9b) 

involve  the  chordwise  forces  described  above  and  must  be  calculated 

separately.  These  preliminary  calculations  require  the  pressure 

coefficients  a and  sectional  generalized  forces  0 obtained 
V(j  -ij 

from  Rowe  program  output.  Program  output  of  generalized  forces  is  of 
the  form 


7-30)  = / /b  ACpj  (x,y)Hi(x,y)  dx  dy 


where  H^(x,y)  is  the  ith  modal  displacement  and  ACp.(x,y)  the 
pressure  distribution  per  unit  q^.  . 

Rederivation  of  (4-9)  with  three-dimensional  subsonic  steadv  air- 
loads produces  the  following  nonlinear  system: 


MPi  ^w  MPi  n n ^v 

(7-31a)  H'n;  -ipS  -ft  - (T-l)  -jp2  [ E E H „•  -J 

J p=l  V=1  J 


n n n 


i*i  Ru^;u  <;v  tt*  - A tq,.« + 
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_ O 

MPi  n n V 

( 7- 3 lb ) n'»]  V - <T-»  TF2  I r,  £,HVU1  i;  -r 


xMPi 

a 


V-l  y*l 


rV 


n n n 


(1  < j < n) 


+ E E E R . q!  ql  -r1  1 

v*l  U=1  i-1  iUV  b 


- MZ  \fo  (/-bACpa  dx]  ft  (y)fv/y)dy  q;  a 
v>=l  v .1 

‘2M  1 1 (x’y)dx]f0  (y)fv  (y)dy  qJ 

u= 1 v— l ry  Tv  j y Yv 

+ 4/Tm7  {/Z(Fv)2f  dv  + 2 E / F"  F®  f dy  a q? 


0 a v 


, o a ® v, 
v=l  Vv  j 


+ E E fl  F"  FT  f dy  q“  q?  } = 0 
y=l  v«l°  K vj  K 


(731c) 


Mi  MPi  n n n 

yi2^)2 -^qJ/CT-D-D^t J,  .1  AW 


’w  \ 


v-l  y-1  i=l  uvl1  b ~b  A 


n n n 


_ O O 

q q 
v 


^ ° _ o 

q q 

v V 

V Mi 


- 1 E E 


v=l  y=l  1=1 


yvjj 


b b i y=l  V-l 
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In  (7-31)  the  notation  0*  . refers  to  the  steady  generalized  forces 

*■  * J 

conputed  directly  by  the  Rowe  program.  Here  2n+l  structural  modes 
must  be  used  to  compute  generalized  forces:  n in  vertical  bending,  n 
in  elastic  torsion,  and  one  rigid  pitching  mode. 

Computation  of  the  chordwise  load  terms  in  ( 7—  31b ) quickly  becomes 
an  unwieldy  task  as  n is  increased,  because  of  the  profusion  of 
numerical  integrations  that  become  necessary  (for  n - 3 there  would 
be  75  integrations).  In  order  to  avoid  extensive  labor  but  still  retain 
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the  most  significant  effects  of  the  potential  chordwise  forces,  the 


following  simplifying  assumptions  are  adopted.  First,  integrals 
appearing  in  the  first  (j  = l)  chordwise  equation  of  ( 7—  31b ) only  are 


retained.  Second,  only  the  pressure  distribution  and  slope  of  the  first 


of  the  elastic  torsion  modes  are  kept;  terms  containing  ACp°  or  f, 

% 

for  u > 1 are  dropped.  Similarly,  all  terms  containing  F"  for 


U 


M > 1 are  removed.  This  approach  should  preserve  the  first-order 
effect  of  steady  chordwise  potential  forces,  yet  only  five  numerical 
spanwise  integrations  will  be  required. 

Actual  computation  of  the  remaining  integrals 


Sq  [/bb  ACp°(x,y)dx]f^  (y)fv  (y)dy 

^ C^bACpJ1(x*y)dxlf(|)1(y)fv1(y)dy 


is  done  by  direct  calculation  of  the  integrals  over  x as  sectional 
generalized  forces  with  the  Rowe  program  at  eleven  spanwise  stations. 
The  spanwise  integration  is  then  carried  out  numerically.  For  the 
three  integrals 

fo  [Fa(y)l2fV;L(y)dy 

fo  Fa(y)F^1<:y)fv1(y)dy 

/0Vj(y)]2fv  (y)dy 

l i 


the  summations  implied  by  (7-18)  are  first  made.  This  step  is  followed 
by  spanwise  numerical  integration. 
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The  nonlinear  solution  scheme  described  in  Appendix  C was  straight- 


forwardly adapted  to  compute  steady  deflections  in  subsonic  three- 
dimensional  flow  from  (7-31).  Since  spanwise  induction  introduces 
coupling  among  the  torsional  modal  equations,  solution  for  linear  dis- 
placements as  an  initial  estimate  becomes  an  nxn  linear  matrix  problem. 

Equation  (4-12)  remains  as  the  3n  x 3n  linear  unsteady  modal 
system  for  stability  about  the  steady  equilibrium  position.  It  now 
requires  the  generalized  force  matrix  to  be  expressed  for  three-dimen- 
sional unsteady  compressible  flow.  As  for  the  steady  case,  the 
generalized  forces  relating  pressures  and  displacements  in  the  vertical 
bending  and  torsion  modes  can  be  computed  directly  by  the  Rowe  program, 
with  reduced  frequency  and  Mach  number  specified.  Direct  insertion  into 
(4-l2a)  and  (4-12c)  is  accomplished  by  the  simple  substitution,  for 
generalized  forces  Q . and  (1  . relating  the  same  two  modes, 

. I f * .If* 


(7-32) 


W -1,1 


Even  more  so  than  for  the  steady  airloads,  complete  inclusion  of 
all  linear  unsteady  potential  chordwise  terms  introduces  a profusion 
of  Integrals.  Formally,  the  linear  unsteady  terms  which  appear  in 
(7-28)  and  (7-23)  enter  into  generalized  forces,  computation  of  which 
involves  spanwise  integrations  of  the  terms'  products  with  fv 
Practically,  simplifying  assumptions  of  the  type  made  for  the  stendv 
chordwise  terms  have  been  made  to  keep  the  number  of  numerical  Integra- 
tions at  a manageable  level.  Accordingly,  only  the  first  chordwise 
equation  of  (4-12b)  will  retain  the  chordwise  force  terms.  Furthermore, 
only  the  chordwise  force  integrals  containing  pressure  contributions 
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and  displacements  for  the  q , q , and  q degrees  of  freedom  will 


‘♦l  W1 


w- 


be  kept.  (Justification  of  this  approximation  is  based  on  the  flutter 
behavior  observed  in  Chapter  VI,  which  revealed  little  participation  by 
the  remaining  vertical  bending  and  torsion  modes.)  The  remaining  linear 
unsteady  terms  from  (7-23)  and  (7-28)  are 

qw  (t) 

(7-33)  Dp(y;t)  = hpV2{f\  ACp^  f^  dx  qj 

qw  (t) 

+ ;-bb  ACpw2f^dx  qJ1  + ;-bb  ACp<01  f<j)1  dx  < q-  (t) 


*♦1  ‘*1 


+ X-bb  Acpa  \ dX  “ + ;-bbACp^  q0,  q*,(t)1 

qw  (t)  qw.  (t) 

+ Snpv2/!^  {F“  F a -r1  + F®  F a 

a a w^  b a w,,  b 


+ K ^ « <u  <o  + K K.  q 


qw  (t) 
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a <J>!  \ 


1 W1  *1  b 


q (t) 


+ F'  FJ  q»  + F®  F.  q ? q.  (t)  > 

<t>2  2 ^1  b ^1  M ^1  ^1 


The  generalized  force  matrix  terms  for  inclusion  in  (4-1 2b)  are  then 
found  to  be 


(7-34) 


Ql+n,l  " M /o'l^bbACpw1(x-y)dx,f«t>1(y)fv1(y)dy  qJT 


f-77  /T-mT  {/f  F®  (y) F (y ) f (v)dvot  +/?  F®  F f dv  q°  F 
k fc  a 0 a Wj  Vj  0 Wj  v^ 
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All  remaining  Q,.  , terms  in  (4-12b)  will  be  zero, 

j +n , i 

For  a given  reduced  frequency,  computation  of  these  generalized 
forces  involves  program  output  for  both  oscillating  and  steady  flow 
conditions.  Sectional  generalized  forces  and  pressure  series  coef- 
ficients output  by  the  Rowe  program  are  used,  in  the  same  manner 
described  above  for  steady  chordwise  loads,  to  compute  (7-34).  Integrals 
for  oscillating  flow,  of  course,  are  complex.  Nine  complex  spanwise 
numerical  integrations  are  needed  for  each  k , whereas  126  would  have 
been  required  without  the  simplifying  approximations.  Two  real  integra- 
tions in  Q,  , involving  steady  pressures  also  appear  in  the 

1+n , I+zn 

steady  displacement  solution. 


B.  Flutter  Calculation  Procedure  and  Results 

Inclusion  of  subsonic  three-dimensional  (3-1))  airloads  eliminates 
a serious  flaw  of  the  strip  theory  loads  used  to  obtain  all  Chapter  VI 
results.  This  was  the  approximate  spanwise  load  distribution,  which  is 
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most  inaccurate  near  the  tip.  Leading-edge  suction  has  been  included 

in  the  same  manner  as  the  Chapter  V,  Section  D,  analysis  of  2-P  flow. 

An  effect  which  was  neglected  in  the  strip  theory  case,  that  of  the 

x-component  of  the  resultant  pressure  force  normal  to  the  deformed 

chord  (cf.  Eqs.  (5-18)  to  (5-20)  is  now  retained.  It  is  accounted  for 

by  the  chordwise  terms  which  are  computed  with  sectional  generalized 

forces.  The  influence  of  induced  drag  upon  both  steady  deformations 

and  flutter  stability  should  now  be  implicitly  included  by  the  modeling 

of  chordwise  loads.  Two  parameters  which  must  be  specified  in  addition 

to  those  mentioned  in  Chapter  IV,  in  order  to  define  a specific  wing 

£ 

are  the  aspect  ratio  (-^)  and  Mach  number. 

The  Chapter  IV  stability  calculation  method  has  been  modified  to 
accept  the  subsonic  3-D  steady  and  unsteady  airloads  derived  in  Section 
A.  Since  these  airloads  are  now  externally  generated,  iteration  of  the 
reduced  frequency  to  find  neutrally  damped  eigenvalues  is  no  longer 
feasible.  Generalized  forces  have  to  be  computed  beforehand,  both  for 
steady  flow  and  for  oscillatory  flow  at  preselected  reduced  frequencies. 
Computation  of  steady  displacements  for  a given  angle  of  attack  and 
flight  speed  is  accomplished  as  before,  now  based  on  Equations  (7-31). 
Then  for  each  of  the  preselected  reduced  frequencies,  the  previously 
calculated  generalized  forces  (7-32)  and  (7-34)  are  input,  the  eigen- 
value determinant  (4-19)  is  assembled,  and  complex  eigenvalues  are 
determined.  For  each  of  these  there  are  an  associated  speed  and  damping 
as  shown  by  (2-25)  and  (2-27).  The  speed  and  damping  of  each  flutter 
mode  are  then  plotted  for  all  k , and  a neutral lv-damped  speed  U_ 

V 

is  determined  by  graphical  interpolation.  The  procedure  is  repeated 


with  a newly  estimated  U ; with  care,  U and  U_  can  be  matched 

e e F 

with  sufficient  accuracy  after  two  such  steps. 

The  expense  in  computer  time  to  execute  the  Rowe  program  and  the 
additional  effort  required  to  prepare  the  chordwise  loads  make  it 
desireable  to  use  as  few  as  possible  k values.  Fortunately,  all  U 
vs.  g interpolation  graphs  proved  to  be  quite  smooth.  Flutter  speeds 
accurate  to  three  significant  digits  (sufficient)  were  reliably  obtained 
The  moderate-aspect-ratio  wing  whose  stability  boundary  for  2-D 
airloads  appears  in  Fig.  6-7  has  been  re-analyzed  using  3-D  incompres- 
sible aerodynamic  theory.  Figure  7-1  shows  the  results,  compared  with 
the  2-D  flutter  calculations.  The  curve  marked  "100%  suction"  was 
computed  with  the  complete  system  (4-12),  (7-32),  and  (7-34)  whereas 
that  marked  "0%  suction"  was  found  by  repeating  the  analysis  after 
removal  of  all  terms  containing  the  singularity  strength  parameters 
F etc.  in  (7-34).  This  latter  result  thus  represents  the  effects  of 


1 


forces  normal  to  the  deflected  airfoil  chord  only  in  the  dynamic  equa- 
tions. The  coupling  of  fore-and-aft  bending  motions  into  flutter  bv 
the  leading  edge  suction  forces  is  therefore  absent  in  the  "0%  suction" 


case. 

The  2-D  stability  boundary  involves  the  aerodvnamic  loads  of 
(4-18),  which  actually  do  not  account  for  suction.  The  effect  of  model- 
ing the  suction  force,  as  derived  in  Section  C of  Chapter  V,  is  shown 
by  Figs.  6-35  and  6-36.  These  plots  suggest  that,  if  a stability 
boundary  had  been  determined  with  suction  accounted  for,  the  2-D  curve 
in  Fig.  7-1  would  show  less  influence  of  steady  deformation  and  would 
not  drop  below  U = 7 . The  solid,  3-D  curve  does,  of  course,  account 
for  suction. 
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Comparison  between  the  2-D  and  3-D  flutter  speeds  is  also  made 
more  difficult  because  they  involve  different  steady  aerodynamics.  In 
order  to  minimize  differences  of  this  kind,  the  steady  tip  deflection 
was  judged  to  be  the  best  common  measure  of  steady  airload  effects. 

Inspection  of  Fig.  7-1  affirms  that  the  influence  of  steadv 
deformations  upon  flutter  is  not  just  governed  by  elastic  bending-torsion 
coupling  but  is  also  sensitive  to  the  manner  in  which  the  potential 
aerodynamic  loads  are  applied  to  the  fore-and-aft  degree  of  freedom. 

Even  though  the  chordwise  force  components  represent  tilting  of  a 
relatively  large,  approximately  vertical  resultant  force  vector,  they 
are  seen  to  have  a significant  stabilizing  influence.  This  conclusion 
follows  from  comparing  the  "100%  suction"  and  "0%  suction"  curves.  Any 
analysis  of  this  type  will  obviously  be  sensitive  to  the  way  in  which 
chordwise  forces  are  accounted  for. 

The  influence  of  compressibility  is  next  explored  by  repeating  the 
foregoing  calculations  with  airloads  computed  for  Mach  numbers  of  .6 
and  .8  . Results  are  shown  in  Fig.  7-2.  Strictly  speaking,  this 
procedure  involves  an  inconsistency,  since  Mach  number  is  held  fixed 
while  velocity  is  freely  varied.  The  type  of  calculations  required  to 
model  properly  a wing  at  high  subsonic  speeds  would  require  iterative 
matching  of  Mach  number  and  speed  Up  for  a given  flight  altitude. 

This  refinement  is  deemed  to  be  excessively  costly  in  both  computer 
time  and  effort.  Nevertheless,  the  results  of  Fig.  7-2  provide  inter- 
esting qualitative  information  on  Mach  number  effects. 

As  M^  is  increased,  the  decrease  in  flutter  speed  with  steady 
deformation  for  the  "0%  suction"  case  becomes  less  pronounced.  The 
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results  for  "100%  suction"  show  less  sensitivity  to  Mach  number.  In 
either  case,  no  adverse  effects  due  to  compressibility  upon  flutter  of 


a lifting  wing  are  revealed,  other  than  those  associated  with  a modest 
decrease  in  with  increasing  . 

When  extending  this  analysis  to  higher  subsonic  Mach  numbers,  it 
must  be  remembered  that  the  only  aerodynamic  loads  being  included  are 
those  which  arise  from  inviscid,  first-order,  small  perturbation  steady 
and  unsteady  theory  for  planar  lifting  surfaces.  Induced  drag  is  present, 
as  discussed  above,  but  all  chordwise  forces  which  arise  from  either 
direct  viscous  shears  or  from  modifications  to  the  pressure  distribution 
due  to  the  presence  of  a boundary  layer  are  not  modeled.  Yet  viscous 
effects  of  this  type  become  increasingly  more  pronounced  as  flight 
speeds  approach  the  transonic  regime  and/or  as  mean  angle  of  attack  is 
increased . 

The  chordwise  laods  associated  with  direct  viscous  shear  should  not 
contain  a significant  unsteady  component,  at  least  at  the  lower  reduced 
frequencies  encountered  here.  It  is  expected  that  they  will  produce  a 
steady  drag  force  aligned  with  the  airfoil  section.  By  contrast,  the 
chordwise  potential  loads  represent  the  horizontal  component  of  a 
relatively  large  resultant  circulatory  lift  vector.  Tilting  of  this 
vector  can  introduce  unsteady  chordwise  loads  of  considerable  magnitude. 
Indeed  their  importance  in  aeroelastic  stability  is  suggested  by 
comparison  of  the  100%  and  0%  suction  curves  in  the  figures. 

Since  viscous  shear  should  cause  a predominantly  steady  drag  loading, 
qualitative  information  regarding  its  effect  on  flutter  can  be  inferred 
from  the  strip  theory  studies  of  Chapter  VI  involving  the  drag  parameter 


C.  For  example,  no  substantial  alteration  of  the  dependence  of  flutter 


speed  upon  steady  lift  should  be  expected  from  including  this  additional 
drag  term.  The  effect  of  unsteady  viscous  contributions  to  aeroelastic 
stability  has  not  been  considered  anywhere  in  this  investigation.  The 
mere  prediction  of  unsteady  viscous  chordwise  loads  is  still  regarded 
as  an  open  question  for  experimental  and  analytical  research,  especiallv 
when  turbulent  boundary  layers  are  involved. 

As  a final  case,  the  large-aspect-ratio  example  of  Fig.  6-13  has 
been  reanalyzed.  Aspect  ratio  is  fixed  at  20  and  Mach  number  at  0. 
Results  are  shown  in  Fig.  7-3.  Interestingly,  the  zero-lift  flutter 
speeds  for  2-D  and  3-D  flows  are  nearly  the  same:  of  course,  strip 
theory  is  expected  to  be  more  accurate  for  the  larger  aspect  ratio. 

The  difference  between  the  2-D  and  3-D  (with  suction)  stability 
boundaries  is  actually  deceiving,  since  (as  already  mentioned)  the  2-D 
results  do  not  contain  the  improved  modeling  of  the  suction  force  from 
Chapter  V.  Figures  (6-37)  and  (6-38)  suggest  the  stabilizing  effect 
which  introducing  suction  would  have  upon  the  2-D  curve  in  Fig.  7-3. 

The  key  observation  here  is  that  the  same  type  of  instability,  involving 
substantial  participation  of  the  fore-and-aft  bending  degree  of  freedom 
observed  for  2-D  aerodynamic  loads,  is  still  observed  after  3-D  aero- 
dynamics are  introduced.  Also  the  influence  of  steady  deformation  upon 
flutter  speed  remains  appreciable.  Removal  of  the  unsteady  suction 
terms  decreases  stability,  much  as  was  observed  for  ^ = 10  in  Fig.  7-1. 
Points  on  the  stability  boundary  for  a < .005  could  not  be  found 
reliably  because  of  the  extremely  light  damping  in  this  region: 
approximate  curves  are  shown  by  dashed  lines. 
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In  conclusion,  the  results  found  In  this  chapter  seem  to  confirm 
that  the  phenomena  analyzed  extensively  in  Chapter  VI  are  still  observed 
when  the  strip  theory  airloads  are  replaced  by  airloads  from  3-D 
aerodynamic  theory.  Furthermore,  the  role  of  chordwise  forces  due  to 
leading  edge  suction  has  been  found  to  increase  the  aeroelastic  stability 
of  a wing  undergoing  steady  deformation  due  to  lift. 
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2 D AIRLOADS 


FIGURE  7-3  Flutter  Speeds  as  Influenced  by  Steady  Deformation;  Large-Aspect-Ratio  Example 


Chapter  VIII 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  following  conclusions  are  drawn  from  the  variety  of  results 

obtained  in  this  investigation. 

1)  The  influence  of  steady  drag  on  flutter  speed  changes  from  favorable 
to  unfavorable  as  aspect  ratio  is  increased.  The  frequency  of  the 
second  transverse  bending  mode  decreases,  tending  toward  the  funda- 
mental bending-torsion  flutter  frequency,  as  aspect  ratio  is 
increased . 

2)  The  prediction  of  the  influence  of  steady  drag  upon  flutter  is  not 
substantially  altered  when  steady  deformations  due  to  lift  are 
considered.  The  major  effect  of  steady  drag  is  to  reduce  divergence 
speed,  especially  for  large  aspect  ratios. 

3)  When  a wing  has  such  a large  aspect  ratio  that  its  fundamental 
fore-and-aft  bending  frequency  is  less  than  the  frequency  of  bending- 
torsion  flutter  at  zero  steady  lift,  an  instability  associated  with 
chordwise  bending  occurs.  The  critical  speeds  are  lower  than  the 
zero-lift  flutter  speed  when  any  steady  lifting  deformation  is 
present.  This  type  of  flutter  can  disappear  at  small  steady  deflec- 
tions when  realistic  structural  damping  is  introduced,  but  for 
reasonably  large  displacements  it  can  still  occur. 

4)  Steady  deformations  decrease  flutter  speed  and  flutter  frequency. 

The  effect  is  most  pronounced  when  the  fundamental  fore-and-aft 
bending  frequency  is  near  the  zero-lift  bending-torsion  flutter 
frequency. 
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5)  The  aeroelastic  phenomena  predicted  using  incompressible  strip- 
theory  airloads  are  also  observed  when  three-dimensional,  compres- 
sible subsonic  airloads  are  employed. 

6)  The  inclusion  of  unsteady  leading-edge  suction  forces  moderates 
the  predicted  decrease  in  flutter  speed  due  to  steady  lifting 
deformat  ions . 

7)  In  subsonic  compressible  flow,  with  unsteady  potential  chordwise 
forces,  there  seem  to  be  only  slight  adverse  effects  of  M on 

SL 

flutter  when  steady  deformations  are  present. 

One  of  the  predictions  of  this  study  is  that  a high  performance 
sailplane  undergoing  a limit  load  factor  pullout  from  a dive  could 
encounter  flutter  or  divergence  at  lower  speeds  than  might  be  antici- 
pated from  a conventional  aeroelastic  analysis.  The  increased  steady 
drag  which  accompanies  higher  C would  reduce  the  divergence  speed 
considerably,  while  the  deformation  of  the  flexible  high-aspect-ratio 
wings  would  change  the  dynamic  aeroelastic  stability  as  well. 

The  present  analysis  could  be  refined  still  further.  Nonlinear 
aerodynamic  effects  deserve  further  attention.  Higher  mean  angles  of 
attack  would  lead  to  increased  importance  of  the  turbulent  boundary 
layer,  culminating  in  separation  (stall),  which  will  alter  the  stabiliz- 
ing contribution  found  due  to  leading-edge  suction  in  attached  flow. 

Vehicles  intended  to  operate  at  the  higher  subsonic  or  low  tran- 
sonic speeds  (certain  RPV's  or  missiles  could  still  be  designed  with 
the  straight  wings  considered  here)  are  expected  to  encounter  various 
phenomena  which  could  greatly  affect  the  chordwise  loads.  Viscous 
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shear,  boundary  layer  modification  of  the  flow,  and  thickness  effects 
with  the  appearance  of  shocks  will  all  modify  the  aeroelastic  behavior. 


1 


Finally,  the  transonic  and  supersonic  flow  regimes,  where  drag 
loads  are  considerably  larger  than  at  subcritical  speeds,  remain  largely 
uninvestigated.  In  these  ranges  different  structural  and  aerodynamic 
configurations  are  likely  to  be  involved. 
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Appendix  A 


THE  ASSUMED  MODES 


1.  Bending 

The  natural  mode  shapes  and  frequencies  in  bending  of  a uniform 
beam  of  length  SL  are  found  by  seeking  homogeneous  simple  harmonic 
solutions  of 


El  w""  + mw  = 0 
x 


Letting 


w = W e 


icot 


the  general  solution  is 


W - A^  sin  J ^ y + A^  cos  J ^ y + A^  sinh  y + A ^ cosh 


where 


o El 
2 x 


Application  of  the  bending  boundary  conditions 


W(0)  = W'(0)  = W "(«,)  = W'"  («,)  = 0 


results  in  a trancendental  equation  for  the  natural  frequency  eigenvalues 


cos 


A 


cosh  — £ 
a 


which  are  for  vertical  bending 
<A-i)  ^ - n!»J  £ J 5 


i - 1,2,3,, 
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PHKCSW.NO  P tat 


The  trancendental  numbers 


for  1 < i < 5 are 


(A- 2) 


Nx  = 0.596864162695 

N2  = 1.494175614274 

N3  = 2.500246946168 

N.  = 3.499989319849 
4 

N5  = 4.500000461516 


and  the  corresponding  eignevectors  yield  the  vertical  bending  natural 
modes 

f sin  N II  - sinh  N.TU 

— , — — ~r: — (sinh  JIN  y - sin  riN.y) 

^osh  TTNi  + cos  IIN^  j 1 

+ cosh  IIN^y  - cos  ITN^y] 


(A- 3) 


expressed  in  orthonormal  form  so  that 

f0  y)dy  = i 

The  modal  property  in  the  steady  equations  (4-9)  is  related  to 

the  modal  integral 


(A-4) 


1 

'o  f 


w 


dy 


"sin  irN^  - sinh  * 

cosh  rnN.  + cos  7rN, 

L i i 


Fore-and-aft  bending  natural  modes  f are  the  same  as  in  (A-3), 

vi 

however  the  natural  frequencies  are 


(A- 5) 
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The  f are  illustrated  in  Fig.  A-l  for  1 < i < 5 . 
w — — 

i 

2.  Torsion 

Natural  mode  shapes  and  frequencies  in  torsion  are  simply 
(A-6)  f^  = sin  H(j-Ss)y 

1 

(A-7)  ^ = IKj-J*) 

which  result  from  the  elementary  Sturm-Liouville  problem  for  torsion 
of  a uniform  rod.  These  modes  are  not  normalized  since 

fo  f^(y)dy  = 35 
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FIGURE  A-l  Assumed  Bending-Mode  Shapes 


Appendix  B 


CALCULATION  OF  MODIFIED  BESSEL  FUNCTIONS 


The  modified  Bessel  functions  K^(s)  and  K^(s)  appearing  In 
the  expression  (5-1)  for  the  generalized  Theodorsen  function  can  be 
computed  by  using  the  following  ascending  power  series  expansions, 
drawn  from  Abramowitz  and  Stegun  (Ref.  26,  Equations  9.6.10,  9.6.11, 
and  9.6.13). 


(B-o  i v<;>  - («)v  ^ 


(B-2) 


K (s)  = - {£n(^s)  + y }l  (s) 
o e o 


(B-3) 


i '4s2  i mM  (ta2)2  , n ,t  . K (hz2)3  . 

+ QTp-  + (2!)'z  ' + (™(jijT'  + 

Kx(s)  = ^(^s)  1 + £n(2 s) Ij(s) 


(hs  2 ) j 


h(hl)  z {<KJ  + l)  + <K1  + 2)}  --TfttnT 

j=o  ' J 


where 


T(v  + 1)  = v! 
4»(1)  = - Y„ 


<p(v)  = - y + Z j , V > 2 
6 .1=1 


Y = 0.5772156649...  (euler’s  constant) 
e 


No  convergence  difficulties  with  the  power  series  expansions  were 
encountered  over  the  range  of  s that  occurred  in  this  investigation. 
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Simple  harmonic  airloads  involving  the  Theordorsen  function  of  reduced 
frequency  C(k)  were  computed  by  the  same  procedure,  simply  by  the 
s = 0 + ik  . 


substitution 


Appendix  C 


SOLUTION  FOR  NONLINEAR  STEADY  DISPLACEMENTS 


With  the  vector  of  generalized  displacements  defined  as 


(C-l)  q ° = {q°  q,°  ,q°,,q°  ,q“  ,..,q°  > 

~ \ Wn  Vn  h \ 


the  nonlinear  equations  (4-9)  can  be  expressed. 


(C-2)  F(q°)  = 0 


where 


(C-3)  F = {F  , ...F  ,F  ,...,F  ,F,  , ...,F.  } 

~~  w ’ w v,  v d> 

1 n 1 n 1 n 


Let  an  initial  estimate  q°,  . be  found  by  solving  the  linear  steady 

w^O) 

equations  (4-20),  then  linearize  F about  q?  by  first-order  Tavlor 

~v  ^ (o; 

series. 


(C 


■4)  _F(£)  - F(,Jo))  + lW(o))l<q‘  - q;o))  + "-O.T. 


The  Jacobian  matrix  J contains  partial  derivatives  of  the  F's  with 
respect  to  the  q°'s  evaluated  at  5 its  elements  are  shown  in 

Fig.  C-l.  Equation  (C-4)  can  be  used  to  solve  for  (q°  . - 9/^)  by 
the  linear  approximation 

<c-5)  J<a (o)> + 'S.W  - 0 

and  the  first  iteration  solution  is 


(C-6) 


i«)  ■ <a<i)  - ’(o)J  +j(o) 
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This  process  is  repeated  until  satisfactory  convergence  is  achieved, 
with  the  general  jth  step  given  by 


(c-7) 

i.U+1)  ■ (2(j+i)  ~~(i>7  * 

In  practice,  more  than  four  or  five  iterations  rarely  were  required, 
for  cases  near  static  divergence  speeds. 
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9F 

w,  MPi  MPi  n n 

Jj.i  = ?q^  = blF-  6ij  + (T_1)  bU7-  Rijyv  <>;u  q£ 


dFw  MPi  n 

Tj,i+n  = 9qJ7  = " vlX  *W  qJv 


3F 

w MPi  n 

7 = 2.  — _ Cr— i i ® r v u P 

j , i+2n  3q°  (T  1J 


r V-i ■ 2 j,  v2!  w \ t>  - »« 


o V, 


9F 


v,  MPia  n 


3+n.i  ^ 


bU 


^ (t"i)  * Hvij  q; 

V=1  J 


3F 


MPi 

Jj«,x+  n ' hF2  sij  - <T-»  rrn£  r E « -*  -• 


MPi_  n n 

bU*  - “ Rijyvq<ji  q4> 

v=l  y=l 


9F 


MPi  n % 


j+n,i+2n  3q 


9Fa 


r1  - - (t-i)  f e»  l7-u+  2 z 


n n 


v,>jb  • - u:x  v.Ex^vi + 2a,ij 


<}>.  MPi  n n 

Jj+2n,i+n  = 9q°  = (T_1)  mF~  2 Z Z R- 

Vi 


V,  n <C 

_ H.  a°  _ y u y.i 

ti^i  b % vi  -ivi  » 


9F, 


Mi 


Jj+2n,i+2n  = 3^f  = - A)  6 

<Pi  J 


MPi  n n „ 

+ (T-l)  -stS  [ E E R — - 


« ° - o 

qv  qv 

-EE  R r~“  r— ~J 


X.  V.  n n 


U lv=ly=l  ^ b b AAW-b 


ij 


1 (i-j) 

(o  (i*J) 

(i  < i < n)  , (i  < j < n) 

FIGURE  C-l  Elements  of  the  Jacobian  Matrix 
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